matlab数据序列的几种滤波器

 

一维数据序列滤波的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('小波滤波后信号 ');

 

 

posted @ 2016-08-05 00:19  MnsterLu  阅读(22948)  评论(2)  编辑  收藏  举报