Matlab: 白噪聲與曲線擬合

在信號處理中常常需要用到曲線擬合，這里介紹一下利用最小二乘擬合一般曲線的方法，並對濾掉信號中白噪聲的方法作些介紹。
為了測試擬合算法的好壞，先模擬出一個信號作為檢驗算法的例子：

1. 用白噪聲產生模擬信號：
   對於理論信號y=y(x)，一般可用rand(size(x))和randn(size(x))生成隨即噪聲信號，兩者的區別在於rand生成的噪聲信號都是正值，而randn生成的噪聲信號則是正負跳躍分布的，所以randn作為白噪聲信號，更符合實際情況：

f0=@(c,x)( (x>=0&x<c(1))*0 + (x>=c(1)&x<c(2))*c(3)/(c(2)-c(1)).*(x-c(1)) + (x>=c(2)&x<c(4)).*( (c(5)-c(3))/(c(4)-c(2))*(x-c(2))+c(3) ) + (x>=c(4)&x<c(6))*c(5)/(c(4)-c(6)).*(x-c(6)) + (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));

2. 最小二乘折線擬合
   考慮到需要擬合的函數是個分段的折線函數，需要首先建立含有固定參數的折線函數的數學模型，算法如下圖：
   
   按照這個算法，用matlab搭建的代碼如下：

% try zigzag fitting
f2=@(c,x)( (x>=0&x<c(1))*0 + (x>=c(1)&x<c(2))*c(3)/(c(2)-c(1)).*(x-c(1)) + (x>=c(2)&x<c(4)).*( (c(5)-c(3))/(c(4)-c(2))*(x-c(2))+c(3) ) + (x>=c(4)&x<c(6))*c(5)/(c(4)-c(6)).*(x-c(6)) + (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 on
plot(x_int,y2,'red --','linewidth',2);
legend('sig','zigzag fitting');

真實參數：1，2，1，5，2，6
擬合參數：1.0237，2.06，1.0107，4.9479，2.1101，6.0005
可以看到，擬合的參數多少和真實的參數存在一些差異，但是已經非常接近。

3. 優化：傅立葉變換降噪
   如果要進一步提高擬合的精度，需要設法降低白噪聲的干擾。因為白噪聲是一種寬譜的干擾，所以常用的帶通濾波處理是不可行的，這里可以考慮對信號進行傅立葉變換，濾掉其中強度較弱的白噪聲頻域成分。

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 direction
sig_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 filter
Af_level=0.01;
Af_lim=Af_level*max(sig_fft_real);
i_fd=find(sig_fft_real<Af_lim);
sig_fft_fit=sig_fft_comp;
sig_fft_fit(i_fd)=0;
figure();
plot(f_double,sig_fft_real_adjust);
xlabel('f(Hz)');
ylabel('A(f)');
xlim([f_double(1),f_double(end)]);
hold on
plot(f_double,Af_lim*ones(size(f_double)),'red --','linewidth',1);
legend('spectrum','Af limit');
% reconstruct the signal with filtered spectrum
sig_fit=ifft(sig_fft_fit);
% perform fitting for the A filtered signal
disp('fit c0 after A filter');
c_fit3=nlinfit(x_int,sig_fit,f2,c0)
y3=f2(c_fit3,x_int);
% compare signal and fitted signal
figure();
plot(x_int,sig,'black',x_int,sig_fit,'red');
hold on
plot(x_int,y3,'green --','linewidth',2);
legend('sig','Fourier fit','zigzag fit');

傅立葉降噪后結果如下：

此時算得的擬合系數是： 1.0677，1.8680， 0.9665，5.0140，1.9736，5.9895 這比降噪前的效果稍好了一些，更貼近與真實的折線系數。但是編程的復雜度上升了很多，在對擬合的精度要求不是太高的情況下，可以不用作傅立葉降噪的處理。

4. 補充：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}\)