1. 程式人生 > >轉載-七種濾波方式的matlab實現

轉載-七種濾波方式的matlab實現

建立兩個混合訊號,便於更好測試濾波器效果。同時用七中濾波方法測試。 混合訊號Mix_Signal_1 = 訊號Signal_Original_1+白噪聲。 七中濾波方法測試matlab實現(轉) - 夏衡 - 槿の伊甸 七中濾波方法測試matlab實現(轉) - 夏衡 - 槿の伊甸

混合訊號Mix_Signal_2 = 訊號Signal_Original_2+白噪聲。 七中濾波方法測試matlab實現(轉) - 夏衡 - 槿の伊甸 七中濾波方法測試matlab實現(轉) - 夏衡 - 槿の伊甸

1.巴特沃斯低通濾波器去噪 巴特沃斯濾波器適合用於訊號和噪聲沒有重疊的情況下。下圖是巴特沃斯對兩個訊號的濾波效果。 七中濾波方法測試matlab實現(轉) - 夏衡 - 槿の伊甸

從圖上可以看出巴特沃斯低通濾波器對訊號一的濾波效果還是可以的,主要是因為有效的訊號最高頻率才30Hz,本程式將50Hz以上的訊號全部濾除,通過的頻率成分中仍然是有白噪聲的。 對於訊號二,濾波後的訊號與沒有加噪聲的訊號相比就有失真了,上升沿和下降沿的高頻訊號被濾除了。

2.FIR低通濾波器去噪 情況同巴特沃斯低通濾波器相似。濾波後的效果如下: 七中濾波方法測試matlab實現(轉) - 夏衡 - 槿の伊甸

  1. 移動平均濾波去噪 濾波效果如下: 七中濾波方法測試matlab實現(轉) - 夏衡 - 槿の伊甸

  2. 中值濾波去噪 七中濾波方法測試matlab實現(轉) - 夏衡 - 槿の伊甸

從上圖可以看出,無論是對訊號一還是對訊號二,中值濾波的濾波效果都是很不錯,特備是對於訊號二,上升沿和下降失真比較的小。

  1. 維納濾波去噪

維納濾波器屬於現代濾波器,傳統的濾波器只能濾除訊號和干擾頻帶沒有重疊的情況,當訊號和干擾頻帶有重疊的時候傳統濾波器將無能為力,這時就需要用到現代濾波器,現代濾波器利用訊號和干擾的統計特徵(如自相關函式、功率譜等)匯出一套最佳估值演算法,然後用硬體或軟體予以實現。 維納濾波是以均方誤差最小(LMS(Least MeanSquare)為準則的,它根據過去觀測值和當前觀測值來估計訊號的當前值,因此它的解形式是系統的傳遞函式或單位脈衝響應。 均方誤差為: 七中濾波方法測試matlab實現(轉) - 夏衡 - 槿の伊甸

維納-霍夫(Wiener-Hopf)方程最小均方誤差下的解為: 七中濾波方法測試matlab實現(轉) - 夏衡 - 槿の伊甸

由於理解不深,對於訊號二,沒有什麼濾波效果 七中濾波方法測試matlab實現(轉) - 夏衡 - 槿の伊甸

  1. 自適應濾波去噪 維納濾波器引數是固定的,適合於平穩隨機訊號。卡爾曼濾波器引數是時變的,適合於非平穩隨機訊號。然而,只有在訊號和噪聲的統計特性先驗已知的情況下,這兩種濾波技術才能獲得最優濾波。在實際應用中,常常無法得到訊號和噪聲統計特性的先驗知識。在這種情況下,自適應濾波技術能夠獲得極佳的濾波效能,因而具有很好的應用價值。 自適應濾波的濾波效果如下: 七中濾波方法測試matlab實現(轉) - 夏衡 - 槿の伊甸

本程式是基於LMS演算法的自適應濾波,從上圖可以看出,濾波效果也是很不錯的,特別是對於訊號二,上升沿有失真,下降沿保持還可以,最要的是得到的波形十分的平滑。由此可見自適應濾波極具使用價值。

  1. 小波去噪

首先看一下小波的去噪效果。 七中濾波方法測試matlab實現(轉) - 夏衡 - 槿の伊甸

對於訊號二,小波的去噪效果非常不錯,雖然得到波形不是很平滑,但是上升沿和下降沿保持的非常高,基本可以看到稜角. %**************************************************************************************** % % 建立兩個訊號Mix_Signal_1 和訊號 Mix_Signal_2 % %***************************************************************************************

Fs = 1000; %取樣率 N = 1000; %取樣點數 n = 0:N-1; t = 0:1/Fs:1-1/Fs; %時間序列 Signal_Original_1 =sin(2pi10t)+sin(2pi20t)+sin(2pi30t); Noise_White_1 = [0.3randn(1,500), rand(1,500)]; %前500點高斯分部白噪聲,後500點均勻分佈白噪聲 Mix_Signal_1 = Signal_Original_1 + Noise_White_1; %構造的混合訊號

Signal_Original_2 = [zeros(1,100), 20ones(1,20), -2ones(1,30), 5ones(1,80), -5ones(1,30), 9ones(1,140), -4ones(1,40), 3ones(1,220), 12ones(1,100), 5ones(1,20), 25ones(1,30), 7 ones(1,190)]; Noise_White_2 = 0.5randn(1,1000); %高斯白噪聲 Mix_Signal_2 = Signal_Original_2 + Noise_White_2; %構造的混合訊號

%**************************************************************************************** % % 訊號Mix_Signal_1 和 Mix_Signal_2 分別作巴特沃斯低通濾波。 % %***************************************************************************************

%混合訊號 Mix_Signal_1 巴特沃斯低通濾波 figure(1); Wc=2*50/Fs; %截止頻率 50Hz [b,a]=butter(4,Wc); Signal_Filter=filter(b,a,Mix_Signal_1);

subplot(4,1,1); %Mix_Signal_1 原始訊號 plot(Mix_Signal_1); axis([0,1000,-4,4]); title('原始訊號 ');

subplot(4,1,2); %Mix_Signal_1 低通濾波濾波後訊號 plot(Signal_Filter); axis([0,1000,-4,4]); title(‘巴特沃斯低通濾波後訊號’);

%混合訊號 Mix_Signal_2 巴特沃斯低通濾波 Wc=2*100/Fs; %截止頻率 100Hz [b,a]=butter(4,Wc); Signal_Filter=filter(b,a,Mix_Signal_2);

subplot(4,1,3); %Mix_Signal_2 原始訊號 plot(Mix_Signal_2); axis([0,1000,-10,30]); title('原始訊號 ');

subplot(4,1,4); %Mix_Signal_2 低通濾波濾波後訊號 plot(Signal_Filter); axis([0,1000,-10,30]); title(‘巴特沃斯低通濾波後訊號’);

%**************************************************************************************** % % 訊號Mix_Signal_1 和 Mix_Signal_2 分別作FIR低通濾波。 % %***************************************************************************************

%混合訊號 Mix_Signal_1 FIR低通濾波 figure(2); F = [0:0.05:0.95]; A = [1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0] ; b = firls(20,F,A); Signal_Filter = filter(b,1,Mix_Signal_1);

subplot(4,1,1); %Mix_Signal_1 原始訊號 plot(Mix_Signal_1); axis([0,1000,-4,4]); title('原始訊號 ');

subplot(4,1,2); %Mix_Signal_1 FIR低通濾波濾波後訊號 plot(Signal_Filter); axis([0,1000,-5,5]); title(‘FIR低通濾波後的訊號’);

%混合訊號 Mix_Signal_2 FIR低通濾波 F = [0:0.05:0.95]; A = [1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0] ; b = firls(20,F,A); Signal_Filter = filter(b,1,Mix_Signal_2); subplot(4,1,3); %Mix_Signal_2 原始訊號 plot(Mix_Signal_2); axis([0,1000,-10,30]); title('原始訊號 ');

subplot(4,1,4); %Mix_Signal_2 FIR低通濾波濾波後訊號 plot(Signal_Filter); axis([0,1000,-10,30]); title(‘FIR低通濾波後的訊號’);

%**************************************************************************************** % % 訊號Mix_Signal_1 和 Mix_Signal_2 分別作移動平均濾波 % %***************************************************************************************

%混合訊號 Mix_Signal_1 移動平均濾波 figure(3); b = [1 1 1 1 1 1]/6; Signal_Filter = filter(b,1,Mix_Signal_1);

subplot(4,1,1); %Mix_Signal_1 原始訊號 plot(Mix_Signal_1); axis([0,1000,-4,4]); title('原始訊號 ');

subplot(4,1,2); %Mix_Signal_1 移動平均濾波後訊號 plot(Signal_Filter); axis([0,1000,-4,4]); title(‘移動平均濾波後的訊號’);

%混合訊號 Mix_Signal_2 移動平均濾波 b = [1 1 1 1 1 1]/6; Signal_Filter = filter(b,1,Mix_Signal_2); subplot(4,1,3); %Mix_Signal_2 原始訊號 plot(Mix_Signal_2); axis([0,1000,-10,30]); title('原始訊號 ');

subplot(4,1,4); %Mix_Signal_2 移動平均濾波後訊號 plot(Signal_Filter); axis([0,1000,-10,30]); title(‘移動平均濾波後的訊號’);

%**************************************************************************************** % % 訊號Mix_Signal_1 和 Mix_Signal_2 分別作中值濾波 % %***************************************************************************************

%混合訊號 Mix_Signal_1 中值濾波 figure(4); Signal_Filter=medfilt1(Mix_Signal_1,10);

subplot(4,1,1); %Mix_Signal_1 原始訊號 plot(Mix_Signal_1); axis([0,1000,-5,5]); title('原始訊號 ');

subplot(4,1,2); %Mix_Signal_1 中值濾波後訊號 plot(Signal_Filter); axis([0,1000,-5,5]); title(‘中值濾波後的訊號’);

%混合訊號 Mix_Signal_2 中值濾波 Signal_Filter=medfilt1(Mix_Signal_2,10); subplot(4,1,3); %Mix_Signal_2 原始訊號 plot(Mix_Signal_2); axis([0,1000,-10,30]); title('原始訊號 ');

subplot(4,1,4); %Mix_Signal_2 中值濾波後訊號 plot(Signal_Filter); axis([0,1000,-10,30]); title(‘中值濾波後的訊號’);

%**************************************************************************************** % % 訊號Mix_Signal_1 和 Mix_Signal_2 分別作維納濾波 % %***************************************************************************************

%混合訊號 Mix_Signal_1 維納濾波 figure(5); Rxx=xcorr(Mix_Signal_1,Mix_Signal_1); %得到混合訊號的自相關函式 M=100; %維納濾波器階數 for i=1:M %得到混合訊號的自相關矩陣 for j=1:M rxx(i,j)=Rxx(abs(j-i)+N); end end Rxy=xcorr(Mix_Signal_1,Signal_Original_1); %得到混合訊號和原訊號的互相關函式 for i=1:M rxy(i)=Rxy(i+N-1); end %得到混合訊號和原訊號的互相關向量 h = inv(rxx)*rxy’; %得到所要涉及的wiener濾波器係數 Signal_Filter=filter(h,1, Mix_Signal_1); %將輸入訊號通過維納濾波器

subplot(4,1,1); %Mix_Signal_1 原始訊號 plot(Mix_Signal_1); axis([0,1000,-5,5]); title('原始訊號 ');

subplot(4,1,2); %Mix_Signal_1 維納濾波後訊號 plot(Signal_Filter); axis([0,1000,-5,5]); title(‘維納濾波後的訊號’);

%混合訊號 Mix_Signal_2 維納濾波 Rxx=xcorr(Mix_Signal_2,Mix_Signal_2); %得到混合訊號的自相關函式 M=500; %維納濾波器階數 for i=1:M %得到混合訊號的自相關矩陣 for j=1:M rxx(i,j)=Rxx(abs(j-i)+N); end end Rxy=xcorr(Mix_Signal_2,Signal_Original_2); %得到混合訊號和原訊號的互相關函式 for i=1:M rxy(i)=Rxy(i+N-1); end %得到混合訊號和原訊號的互相關向量 h=inv(rxx)*rxy’; %得到所要涉及的wiener濾波器係數 Signal_Filter=filter(h,1, Mix_Signal_2); %將輸入訊號通過維納濾波器

subplot(4,1,3); %Mix_Signal_2 原始訊號 plot(Mix_Signal_2); axis([0,1000,-10,30]); title('原始訊號 ');

subplot(4,1,4); %Mix_Signal_2 維納濾波後訊號 plot(Signal_Filter); axis([0,1000,-10,30]); title(‘維納濾波後的訊號’);

%**************************************************************************************** % % 訊號Mix_Signal_1 和 Mix_Signal_2 分別作自適應濾波 % %***************************************************************************************

%混合訊號 Mix_Signal_1 自適應濾波 figure(6); N=1000; %輸入訊號抽樣點數N k=100; %時域抽頭LMS演算法濾波器階數 u=0.001; %步長因子

%設定初值 yn_1=zeros(1,N); %output signal yn_1(1:k)=Mix_Signal_1(1:k); %將輸入訊號SignalAddNoise的前k個值作為輸出yn_1的前k個值 w=zeros(1,k); %設定抽頭加權初值 e=zeros(1,N); %誤差訊號

%用LMS演算法迭代濾波 for i=(k+1):N XN=Mix_Signal_1((i-k+1):(i)); yn_1(i)=wXN’; e(i)=Signal_Original_1(i)-yn_1(i); w=w+2u*e(i)*XN; end

subplot(4,1,1); plot(Mix_Signal_1); %Mix_Signal_1 原始訊號 axis([k+1,1000,-4,4]); title(‘原始訊號’);

subplot(4,1,2); plot(yn_1); %Mix_Signal_1 自適應濾波後訊號 axis([k+1,1000,-4,4]); title(‘自適應濾波後訊號’);

%混合訊號 Mix_Signal_2 自適應濾波 N=1000; %輸入訊號抽樣點數N k=500; %時域抽頭LMS演算法濾波器階數 u=0.000011; %步長因子

%設定初值 yn_1=zeros(1,N); %output signal yn_1(1:k)=Mix_Signal_2(1:k); %將輸入訊號SignalAddNoise的前k個值作為輸出yn_1的前k個值 w=zeros(1,k); %設定抽頭加權初值 e=zeros(1,N); %誤差訊號

%用LMS演算法迭代濾波 for i=(k+1):N XN=Mix_Signal_2((i-k+1):(i)); yn_1(i)=wXN’; e(i)=Signal_Original_2(i)-yn_1(i); w=w+2u*e(i)*XN; end

subplot(4,1,3); plot(Mix_Signal_2); %Mix_Signal_1 原始訊號 axis([k+1,1000,-10,30]); title(‘原始訊號’);

subplot(4,1,4); plot(yn_1); %Mix_Signal_1 自適應濾波後訊號 axis([k+1,1000,-10,30]); title(‘自適應濾波後訊號’);

%**************************************************************************************** % % 訊號Mix_Signal_1 和 Mix_Signal_2 分別作小波濾波 % %***************************************************************************************

%混合訊號 Mix_Signal_1 小波濾波 figure(7); subplot(4,1,1); plot(Mix_Signal_1); %Mix_Signal_1 原始訊號 axis([0,1000,-5,5]); title('原始訊號 ');

subplot(4,1,2); [xd,cxd,lxd] = wden(Mix_Signal_1,‘sqtwolog’,‘s’,‘one’,2,‘db3’); plot(xd); %Mix_Signal_1 小波濾波後訊號 axis([0,1000,-5,5]); title('小波濾波後訊號 ');

%混合訊號 Mix_Signal_2 小波濾波 subplot(4,1,3); plot(Mix_Signal_2); %Mix_Signal_2 原始訊號 axis([0,1000,-10,30]); title('原始訊號 ');

subplot(4,1,4); [xd,cxd,lxd] = wden(Mix_Signal_2,‘sqtwolog’,‘h’,‘sln’,3,‘db3’); plot(xd); %Mix_Signal_2 小波濾波後訊號 axis([0,1000,-10,30]); title('小波濾波後訊號 ');