在信号处理中常常需要用到曲线拟合,这里介绍一下利用最小二乘拟合一般曲线的方法,并对滤掉信号中白噪声的方法作些介绍。
为了测试拟合算法的好坏,先模拟出一个信号作为检验算法的例子:- 用白噪声产生模拟信号: 对于理论信号y=y(x),一般可用rand(size(x))和randn(size(x))生成随即噪声信号,两者的区别在于rand生成的噪声信号都是正值,而randn生成的噪声信号则是正负跳跃分布的,所以randn作为白噪声信号,更符合实际情况:
f0=@(c,x)( (x>=0&x=c(1)&x =c(2)&x =c(4)&x =c(6))*0 );disp('real c0');c0=[1, 2, 1, 5, 2, 6]x_int=0:0.002:10;y_int=f0(c0,x_int);%(x_int, y_int) is perfect zigzag signal %sig=y_int+0.5*rand(size(x_int));sig=y_int+0.5*randn(size(x_int));
- 最小二乘折线拟合 考虑到需要拟合的函数是个分段的折线函数,需要首先建立含有固定参数的折线函数的数学模型,算法如下图: 按照这个算法,用matlab搭建的代码如下:
% try zigzag fittingf2=@(c,x)( (x>=0&x真实参数:1,2,1,5,2,6 拟合参数:1.0237,2.06,1.0107,4.9479,2.1101,6.0005 可以看到,拟合的参数多少和真实的参数存在一些差异,但是已经非常接近。=c(1)&x =c(2)&x =c(4)&x =c(6))*0 );c0=[1.1, 1.5, 1.8, 5.4, 2.5, 5.6];c_fit=nlinfit(x_int,sig,f2,c0);y2=f2(c_fit,x_int);figure();plot(x_int,sig,'blue');hold onplot(x_int,y2,'red --','linewidth',2);legend('sig','zigzag fitting');
- 优化:傅立叶变换降噪 如果要进一步提高拟合的精度,需要设法降低白噪声的干扰。因为白噪声是一种宽谱的干扰,所以常用的带通滤波处理是不可行的,这里可以考虑对信号进行傅立叶变换,滤掉其中强度较弱的白噪声频域成分。
Fs=1/(x_int(2)-x_int(1));nfft=length(sig);sig_fft_comp=fft(sig);sig_fft_real=2*abs(sig_fft_comp)/nfft;% adjust the distribution of spectrum according to double frequency directionsig_fft_real_adjust=[sig_fft_real(round(nfft/2+1):end),sig_fft_real(1:round(nfft/2))];f_double=linspace(-Fs/2,Fs/2,nfft);% apply the A(f) strength filterAf_level=0.01;Af_lim=Af_level*max(sig_fft_real);i_fd=find(sig_fft_real
傅立叶降噪后结果如下:
此时算得的拟合系数是: 1.0677,1.8680, 0.9665,5.0140,1.9736,5.9895 这比降噪前的效果稍好了一些,更贴近与真实的折线系数。但是编程的复杂度上升了很多,在对拟合的精度要求不是太高的情况下,可以不用作傅立叶降噪的处理。- 补充:matlab多项式拟合函数(polyfit) [p,s,mu]=polyfit(x,y,n) x,y是被拟合的离散曲线点,n是需要拟合的多项式次数(默认的多项式是幂级数形式的),其中p是个多项式各次项的系数,是按照指数从高到低排列的。mu(1)是y的平均值,mu(2)是单位标准偏差(unit standard deviation,可缩写成STD)\(SDT=\frac{y-mean(y)}{\sigma}\)