一維數據序列濾波的matlab代碼，
其實和之前做的圖像濾波大同小異，
只是圖像的噪聲情況復雜得多，
而且是二維的。
做這個主要是手上有個心電的的mens傳感器，
藍牙把數據傳過來做一個數據的100Hz左右的帶通濾波，
用butterworht做個帶通濾波，
再用c語言重構一下。
1 %**************************************************************************************** 2 % 3 % 創建兩個信號Mix_Signal_1 和信號 Mix_Signal_2 4 % 5 %*************************************************************************************** 6 7 Fs = 1000; %采樣率 8 N = 1000; %采樣點數 9 n = 0:N-1; 10 t = 0:1/Fs:1-1/Fs; %時間序列 11 Signal_Original_1 =sin(2*pi*10*t)+sin(2*pi*20*t)+sin(2*pi*30*t); 12 Noise_White_1 = [0.3*randn(1,500), rand(1,500)]; %前500點高斯分部白噪聲，后500點均勻分布白噪聲 13 Mix_Signal_1 = Signal_Original_1 + Noise_White_1; %構造的混合信號 14 15 Signal_Original_2 = [zeros(1,100), 20*ones(1,20), -2*ones(1,30), 5*ones(1,80), -5*ones(1,30), 9*ones(1,140), -4*ones(1,40), 3*ones(1,220), 12*ones(1,100), 5*ones(1,20), 25*ones(1,30), 7 *ones(1,190)]; 16 Noise_White_2 = 0.5*randn(1,1000); %高斯白噪聲 17 Mix_Signal_2 = Signal_Original_2 + Noise_White_2; %構造的混合信號
一、butterworth filter
1 %**************************************************************************************** 2 % 3 % 信號Mix_Signal_1 和 Mix_Signal_2 分別作巴特沃斯低通濾波。 4 % 5 %*************************************************************************************** 6 7 %混合信號 Mix_Signal_1 巴特沃斯低通濾波 8 figure(1); 9 Wc=2*50/Fs; %截止頻率 50Hz 10 [b,a]=butter(4,Wc); 11 Signal_Filter=filter(b,a,Mix_Signal_1); 12 13 subplot(4,1,1); %Mix_Signal_1 原始信號 14 plot(Mix_Signal_1); 15 axis([0,1000,-4,4]); 16 title('原始信號 '); 17 18 subplot(4,1,2); %Mix_Signal_1 低通濾波濾波后信號 19 plot(Signal_Filter); 20 axis([0,1000,-4,4]); 21 title('巴特沃斯低通濾波后信號'); 22 23 %混合信號 Mix_Signal_2 巴特沃斯低通濾波 24 Wc=2*100/Fs; %截止頻率 100Hz 25 [b,a]=butter(4,Wc); 26 Signal_Filter=filter(b,a,Mix_Signal_2); 27 28 subplot(4,1,3); %Mix_Signal_2 原始信號 29 plot(Mix_Signal_2); 30 axis([0,1000,-10,30]); 31 title('原始信號 '); 32 33 subplot(4,1,4); %Mix_Signal_2 低通濾波濾波后信號 34 plot(Signal_Filter); 35 axis([0,1000,-10,30]); 36 title('巴特沃斯低通濾波后信號');
二、FIR
1 %**************************************************************************************** 2 % 3 % 信號Mix_Signal_1 和 Mix_Signal_2 分別作FIR低通濾波。 4 % 5 %*************************************************************************************** 6 7 %混合信號 Mix_Signal_1 FIR低通濾波 8 figure(2); 9 F = [0:0.05:0.95]; 10 A = [1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0] ; 11 b = firls(20,F,A); 12 Signal_Filter = filter(b,1,Mix_Signal_1); 13 14 subplot(4,1,1); %Mix_Signal_1 原始信號 15 plot(Mix_Signal_1); 16 axis([0,1000,-4,4]); 17 title('原始信號 '); 18 19 subplot(4,1,2); %Mix_Signal_1 FIR低通濾波濾波后信號 20 plot(Signal_Filter); 21 axis([0,1000,-5,5]); 22 title('FIR低通濾波后的信號'); 23 24 %混合信號 Mix_Signal_2 FIR低通濾波 25 F = [0:0.05:0.95]; 26 A = [1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0] ; 27 b = firls(20,F,A); 28 Signal_Filter = filter(b,1,Mix_Signal_2); 29 subplot(4,1,3); %Mix_Signal_2 原始信號 30 plot(Mix_Signal_2); 31 axis([0,1000,-10,30]); 32 title('原始信號 '); 33 34 subplot(4,1,4); %Mix_Signal_2 FIR低通濾波濾波后信號 35 plot(Signal_Filter); 36 axis([0,1000,-10,30]); 37 title('FIR低通濾波后的信號');
三、移動平均濾波
1 %**************************************************************************************** 2 % 3 % 信號Mix_Signal_1 和 Mix_Signal_2 分別作移動平均濾波 4 % 5 %*************************************************************************************** 6 7 %混合信號 Mix_Signal_1 移動平均濾波 8 figure(3); 9 b = [1 1 1 1 1 1]/6; 10 Signal_Filter = filter(b,1,Mix_Signal_1); 11 12 subplot(4,1,1); %Mix_Signal_1 原始信號 13 plot(Mix_Signal_1); 14 axis([0,1000,-4,4]); 15 title('原始信號 '); 16 17 subplot(4,1,2); %Mix_Signal_1 移動平均濾波后信號 18 plot(Signal_Filter); 19 axis([0,1000,-4,4]); 20 title('移動平均濾波后的信號'); 21 22 %混合信號 Mix_Signal_2 移動平均濾波 23 b = [1 1 1 1 1 1]/6; 24 Signal_Filter = filter(b,1,Mix_Signal_2); 25 subplot(4,1,3); %Mix_Signal_2 原始信號 26 plot(Mix_Signal_2); 27 axis([0,1000,-10,30]); 28 title('原始信號 '); 29 30 subplot(4,1,4); %Mix_Signal_2 移動平均濾波后信號 31 plot(Signal_Filter); 32 axis([0,1000,-10,30]); 33 title('移動平均濾波后的信號');
四、中值濾波
1 %**************************************************************************************** 2 % 3 % 信號Mix_Signal_1 和 Mix_Signal_2 分別作中值濾波 4 % 5 %*************************************************************************************** 6 7 %混合信號 Mix_Signal_1 中值濾波 8 figure(4); 9 Signal_Filter=medfilt1(Mix_Signal_1,10); 10 11 subplot(4,1,1); %Mix_Signal_1 原始信號 12 plot(Mix_Signal_1); 13 axis([0,1000,-5,5]); 14 title('原始信號 '); 15 16 subplot(4,1,2); %Mix_Signal_1 中值濾波后信號 17 plot(Signal_Filter); 18 axis([0,1000,-5,5]); 19 title('中值濾波后的信號'); 20 21 %混合信號 Mix_Signal_2 中值濾波 22 Signal_Filter=medfilt1(Mix_Signal_2,10); 23 subplot(4,1,3); %Mix_Signal_2 原始信號 24 plot(Mix_Signal_2); 25 axis([0,1000,-10,30]); 26 title('原始信號 '); 27 28 subplot(4,1,4); %Mix_Signal_2 中值濾波后信號 29 plot(Signal_Filter); 30 axis([0,1000,-10,30]); 31 title('中值濾波后的信號');
五、維納濾波
1 %**************************************************************************************** 2 % 3 % 信號Mix_Signal_1 和 Mix_Signal_2 分別作維納濾波 4 % 5 %*************************************************************************************** 6 7 %混合信號 Mix_Signal_1 維納濾波 8 figure(5); 9 Rxx=xcorr(Mix_Signal_1,Mix_Signal_1); %得到混合信號的自相關函數 10 M=100; %維納濾波器階數 11 for i=1:M %得到混合信號的自相關矩陣 12 for j=1:M 13 rxx(i,j)=Rxx(abs(j-i)+N); 14 end 15 end 16 Rxy=xcorr(Mix_Signal_1,Signal_Original_1); %得到混合信號和原信號的互相關函數 17 for i=1:M 18 rxy(i)=Rxy(i+N-1); 19 end %得到混合信號和原信號的互相關向量 20 h = inv(rxx)*rxy'; %得到所要涉及的wiener濾波器系數 21 Signal_Filter=filter(h,1, Mix_Signal_1); %將輸入信號通過維納濾波器 22 23 subplot(4,1,1); %Mix_Signal_1 原始信號 24 plot(Mix_Signal_1); 25 axis([0,1000,-5,5]); 26 title('原始信號 '); 27 28 subplot(4,1,2); %Mix_Signal_1 維納濾波后信號 29 plot(Signal_Filter); 30 axis([0,1000,-5,5]); 31 title('維納濾波后的信號'); 32 33 %混合信號 Mix_Signal_2 維納濾波 34 Rxx=xcorr(Mix_Signal_2,Mix_Signal_2); %得到混合信號的自相關函數 35 M=500; %維納濾波器階數 36 for i=1:M %得到混合信號的自相關矩陣 37 for j=1:M 38 rxx(i,j)=Rxx(abs(j-i)+N); 39 end 40 end 41 Rxy=xcorr(Mix_Signal_2,Signal_Original_2); %得到混合信號和原信號的互相關函數 42 for i=1:M 43 rxy(i)=Rxy(i+N-1); 44 end %得到混合信號和原信號的互相關向量 45 h=inv(rxx)*rxy'; %得到所要涉及的wiener濾波器系數 46 Signal_Filter=filter(h,1, Mix_Signal_2); %將輸入信號通過維納濾波器 47 48 subplot(4,1,3); %Mix_Signal_2 原始信號 49 plot(Mix_Signal_2); 50 axis([0,1000,-10,30]); 51 title('原始信號 '); 52 53 subplot(4,1,4); %Mix_Signal_2 維納濾波后信號 54 plot(Signal_Filter); 55 axis([0,1000,-10,30]); 56 title('維納濾波后的信號');
六、自適應濾波
1 %**************************************************************************************** 2 % 3 % 信號Mix_Signal_1 和 Mix_Signal_2 分別作自適應濾波 4 % 5 %*************************************************************************************** 6 7 %混合信號 Mix_Signal_1 自適應濾波 8 figure(6); 9 N=1000; %輸入信號抽樣點數N 10 k=100; %時域抽頭LMS算法濾波器階數 11 u=0.001; %步長因子 12 13 %設置初值 14 yn_1=zeros(1,N); %output signal 15 yn_1(1:k)=Mix_Signal_1(1:k); %將輸入信號SignalAddNoise的前k個值作為輸出yn_1的前k個值 16 w=zeros(1,k); %設置抽頭加權初值 17 e=zeros(1,N); %誤差信號 18 19 %用LMS算法迭代濾波 20 for i=(k+1):N 21 XN=Mix_Signal_1((i-k+1):(i)); 22 yn_1(i)=w*XN'; 23 e(i)=Signal_Original_1(i)-yn_1(i); 24 w=w+2*u*e(i)*XN; 25 end 26 27 subplot(4,1,1); 28 plot(Mix_Signal_1); %Mix_Signal_1 原始信號 29 axis([k+1,1000,-4,4]); 30 title('原始信號'); 31 32 subplot(4,1,2); 33 plot(yn_1); %Mix_Signal_1 自適應濾波后信號 34 axis([k+1,1000,-4,4]); 35 title('自適應濾波后信號'); 36 37 %混合信號 Mix_Signal_2 自適應濾波 38 N=1000; %輸入信號抽樣點數N 39 k=500; %時域抽頭LMS算法濾波器階數 40 u=0.000011; %步長因子 41 42 %設置初值 43 yn_1=zeros(1,N); %output signal 44 yn_1(1:k)=Mix_Signal_2(1:k); %將輸入信號SignalAddNoise的前k個值作為輸出yn_1的前k個值 45 w=zeros(1,k); %設置抽頭加權初值 46 e=zeros(1,N); %誤差信號 47 48 %用LMS算法迭代濾波 49 for i=(k+1):N 50 XN=Mix_Signal_2((i-k+1):(i)); 51 yn_1(i)=w*XN'; 52 e(i)=Signal_Original_2(i)-yn_1(i); 53 w=w+2*u*e(i)*XN; 54 end 55 56 subplot(4,1,3); 57 plot(Mix_Signal_2); %Mix_Signal_1 原始信號 58 axis([k+1,1000,-10,30]); 59 title('原始信號'); 60 61 subplot(4,1,4); 62 plot(yn_1); %Mix_Signal_1 自適應濾波后信號 63 axis([k+1,1000,-10,30]); 64 title('自適應濾波后信號');
七、小波濾波
1 %**************************************************************************************** 2 % 3 % 信號Mix_Signal_1 和 Mix_Signal_2 分別作小波濾波 4 % 5 %*************************************************************************************** 6 7 %混合信號 Mix_Signal_1 小波濾波 8 figure(7); 9 subplot(4,1,1); 10 plot(Mix_Signal_1); %Mix_Signal_1 原始信號 11 axis([0,1000,-5,5]); 12 title('原始信號 '); 13 14 subplot(4,1,2); 15 [xd,cxd,lxd] = wden(Mix_Signal_1,'sqtwolog','s','one',2,'db3'); 16 plot(xd); %Mix_Signal_1 小波濾波后信號 17 axis([0,1000,-5,5]); 18 title('小波濾波后信號 '); 19 20 %混合信號 Mix_Signal_2 小波濾波 21 subplot(4,1,3); 22 plot(Mix_Signal_2); %Mix_Signal_2 原始信號 23 axis([0,1000,-10,30]); 24 title('原始信號 '); 25 26 subplot(4,1,4); 27 [xd,cxd,lxd] = wden(Mix_Signal_2,'sqtwolog','h','sln',3,'db3'); 28 plot(xd); %Mix_Signal_2 小波濾波后信號 29 axis([0,1000,-10,30]); 30 title('小波濾波后信號 ');
