振动信号背景噪声滤波算法

@毛毛

drownedfish 发表于 2016-9-9 14:26
这种噪声怎么去掉？没有启动实验台之前，传感器自带的噪声和环境噪声

有这种专门的处理软件吧。。

野菊

@毛毛 发表于 2016-9-10 21:00
有这种专门的处理软件吧。。

什么软件？

猫头鹰先生

王利明 发表于 2016-2-15 06:14
抑制噪声通常以什么为目标函数比较呢，峭度吗，还是熵，还是什么？？

SNR and AE, etc..

猫头鹰先生

OK. Such as Wiener filtering；Butterworth LPF；general chebyshev filtering；chaos pass filtering；median filtering algorithm；wavelet threshold de-noising；Noise Adaptive filtering.

Some filtering examples and their MATLAB codes are as follows:

%****************************************************************************************
%  
%                      创建两个信号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(2*pi*10*t)+sin(2*pi*20*t)+sin(2*pi*30*t);
Noise_White_1    = [0.3*randn(1,500), rand(1,500)];           %前500点高斯分部白噪声，后500点均匀分布白噪声
Mix_Signal_1   = Signal_Original_1 + Noise_White_1;        %构造的混合信号

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)];
Noise_White_2     =  0.5*randn(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)=w*XN';
        e(i)=Signal_Original_1(i)-yn_1(i);
        w=w+2*u*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)=w*XN';
        e(i)=Signal_Original_2(i)-yn_1(i);
        w=w+2*u*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('小波滤波后信号 ');

Pseudo-lover

猫头鹰先生 发表于 2016-9-12 09:40
OK. Such as Wiener filtering；Butterworth LPF；general chebyshev filtering；chaos pass filtering；me ...

太感谢了  刚运行了一下  很棒！！！给力！！！

Pseudo-lover

Noise_White_1    = [0.3*randn(1,500), rand(1,500)];
不都是随机点吗？？？

Pseudo-lover

再补充一下：MATLAB中通过butter函数构造滤波器

[B,A] = BUTTER(N,Wn) designs an Nth order lowpass digital Butterworth filter and returns the filter coefficients in length N+1 vectors B (numerator) and A (denominator).

The cutoff frequency Wn must be 0.0 < Wn < 1.0, with 1.0 corresponding to half the sample rate.

If Wn is a two-element vector, Wn = [W1 W2], BUTTER returns an order 2N bandpass filter with passband W1 < W < W2.
    [B,A] = BUTTER(N,Wn,'high') designs a highpass filter.
    [B,A] = BUTTER(N,Wn,'low') designs a lowpass filter.
    [B,A] = BUTTER(N,Wn,'stop') is a bandstop filter if Wn = [W1 W2].

当设计低通和高通时，Wn是一个值，表示截止频率；当设计带通和带阻时，Wn是一个二个元素的数组，表示通带或阻带的上下截止频率。频率的归一化是对fs/2进行归一。

例如，fs=10000，则fs/2=5000，500归一化为0.1，1000归一化为0.2。

fs=100;%采样频率为100Hz
fc=30;%截止频率为30Hz
[b,a]=butter(4,0.6,'high');
[h,w]=freqz(b,a);
plot(w/pi*fs/2,abs(h)); grid;
title('Amplitude Response');
xlabel('Frequency (Hz)'); ylabel('Amplitude');
figure;
x=randn(1,fs*6);
y=filter(b,a,x);
subplot 211; plot(x); title('Original Signal');
subplot 212; plot(y); title('Output of High-pass filter');
其中信号是用随机数产生的，长6s。笫1张图是高通滤波器的幅值响应曲线，笫2张图是信号滤波前后的波形。

[b,a]=butter(n,Wn)，根据阶数n和截止频率Wn计算ButterWorth滤波器分子分母系数（b为分子系数的矢量形式，a为分母系数的矢量形式）

butter函数是求Butterworth数字滤波器的系数，在求出系数后对信号进行滤波时用filter函数。
说白了，设计滤波器就是设计滤波器系数[B,A]。
[B,A] = BUTTER(N,Wn,'high') ---用来设计高通滤波器
[B,A] = BUTTER(N,Wn,'low') designs a lowpass filter.--低通滤波器
[B,A] = BUTTER(N,Wn)--带通滤波器
N是滤波器的阶数，不熟的话，大概取个整数就可以了。Wn的确定跟你的采样频率Fs有关。
对于原始信号x。
比如说你的采样频率Fs=1000Hz，设计一个8阶、通带为100-200Hz的带通滤波器：
[b,a]=butter(8,[0.2 0.4])=butter(8,[100/(1000/2) 200/(1000/2) ])
这里Fa=Fs/2，Fa是分析频率
得到滤波器系数后，就可以直接用了。
y=filter(B,A,x)

filter(ones(1,5)/5,1,data)
这里b=(1/5 1/5 1/5 1/5 1/5), a=1。
好像有这么个说法：
a=1是FIR----------有限冲激响应滤波器
a不等于1是IIR----无限冲激响应滤波器

Triste

感谢楼上分享的好东西  学习了  谢谢大家

@毛毛

猫头鹰先生 发表于 2016-9-12 09:40
OK. Such as Wiener filtering；Butterworth LPF；general chebyshev filtering；chaos pass filtering；me ...

hi 是否有关于这几种滤波方式的介绍的资料？

sovereign

谢谢大家分享