信号处理(1) –常用信号平滑去噪的方法

前言:最近研究汽车碰撞的加速度信号,在信号的采集过程中难免遇到噪音,导致信号偏差,为了更好的反映系统情况,故常需要信号去噪,本文分享一些
常用信号平滑去噪的方法。

关键字:信号;去噪;Matlab


信号在实际测量中,难免会混入各种噪声。通常我们希望去除高频的随机噪声,或者是偏离正常测量太大的离群误差,以获得低频的测量数据。下面介绍几种常用的信号平滑去噪的方法。


1、移动平均法

滑动平均法(moving average)也叫做移动平均法、平均法、移动平均值滤波法等等,是一种时间域思想上的信号光滑方法。算法思路为,将该点附近的采样点做算数平均,作为这个点光滑后的值。

信号处理(1) --常用信号平滑去噪的方法

一般窗口为对称窗口,防止出现相位偏差。窗口一般为奇数。

以3点平均(窗口长度为3)公式为例,原数据为x,平滑后的数据为y:

y(n)=1/3∗(x(n−1)+x(n)+x(n+1))

对y(n)和y(n+1)相减,可以得到另一种计算形式:

信号处理(1) --常用信号平滑去噪的方法


2、Matlab内自带函数实现移动平均法

matlab有两个函数实现滑动平均法,一个是smoothdata()函数,一个是movmean()函数。

以窗口长度为5为例,smoothdata()函数调用方法为:

y = smoothdata( x , 'movmean' , 5 );

但是这个smoothdata函数实际上是调用了movmean()函数。所以如果直接使用的话,直接用movmean()会更快。

movmean()函数的调用方法为:

y = movmean( x , 5 );

下面以一个加噪声的正弦信号为例:

%移动平均滤波Nber_window = 3;%窗口长度(最好为奇数)t = 0:0.02:3;A = cos(2*pi*2*t)+0.5*rand(size(t));B1 = movmean(A,Nber_window);figure(1)plot(t,A,t,B1)

信号处理(1) --常用信号平滑去噪的方法


3、利用卷积函数conv()实现移动平均法

根据之前的公式,我们可以看到

y(n)=1/3∗(x(n−1)+x(n)+x(n+1))

就相当于一个对x和向量[1/3 1/3 1/3]做卷积。可以验证:

%移动平均滤波Nber_window = 3;%窗口长度(最好为奇数)t = 0:0.02:3;A = cos(2*pi*2*t)+0.5*rand(size(t));B1 = movmean(A,Nber_window);
F_average = 1/N_window*ones(1,N_window);%构建卷积核B2 = conv(A,F_average,'same');%利用卷积的方法计算plot(t,B1,t,B2)

信号处理(1) --常用信号平滑去噪的方法

中间部分两者完全一致,但是两端有所差别。主要是因为,movmean()函数在处理边缘时,采用减小窗口的方式,而conv()相当于在两端补零。所以如何处理边缘也是值得注意的。


4、利用filter滤波函数实现移动平均法

首先介绍一下Z变换。以向前的滑动平均为例(这里中间值不是n而是n+1,所以相位会移动)。

y(n)=1/3∗(x(n)+x(n+1)+x(n+2))

它的Z变换可以简单的理解为,把x(n+k)替换为z^(-k),即

因此对于filter滤波函数,输入的格式为:

y

其中b和a的定义为:

信号处理(1) --常用信号平滑去噪的方法

其中a1必须为1。所以对应的移动平均滤波可以表示为:

y = filter( [1/3,1/3,1/3] , [1] , x );

它和下面代码的是等价的(在边缘上的处理方式有所不同)

y = movmean( x , [2,0] );

代表有偏移的滑动平均滤波,2是向后2个点的意思,0是向前0个点的意思。

因为 filter滤波器使用有偏移的向后滤波。滤波后,相位会发生改变。所以通常采用零相位滤波器进行滤波,matlab内的函数为filtfilt()。原理从函数名字上就可以体现出来,就是先正常滤波,之后再将信号从后向前再次滤波。正滤一遍反滤一遍,使得相位偏移等于0。

信号处理(1) --常用信号平滑去噪的方法


5、移动平均的幅频响应

幅频响应可以通过之前4得到的H(z)函数来得到,在单位圆上采样,也就是把z替换为e^iw。

以中心窗口为例,

信号处理(1) --常用信号平滑去噪的方法

H(iw)的绝对值就是该滤波方法的幅频响应。以3点滤波为例,matlab代码为:

%计算频率响应N_window = 3;w = linspace(0,pi,400);H_iw = zeros(N_window,length(w));n=0;for k = -(N_window-1)/2:1:(N_window-1)/2    n = n+1;    H_iw(n,:) =1/3* exp(w.*1i).^(-k);%公式(e^iw)^(-k)endH_iw_Sum = abs(sum(H_iw,1));%相加figure()plot(w/2/)

信号处理(1) --常用信号平滑去噪的方法

由于H变换在单位圆上的特性相当于傅里叶变换,所以上面代码等效于下面计算方法:

%计算频率响应N_window = 3;Y=zeros(400,1);Y(1:N_window) = 1/3;%设置滤波器Y_F = fft(Y);figure()plot(linspace(0,0.5,200),abs(Y_F(1:200)));

matlab也有自带的函数来看频率特性,freqz(),推荐使用这种。

其中,归一化频率等于信号频率除以采样频率f/Fs,采样频率等于时间采样间隔的倒数1/dt。对比不同窗口长度的幅频响应,可以看到:

1)平均所采用的点数越多,高频信号的滤波效果越好。

2)3点平均对于1/3频率的信号滤波效果最好,5点平均对1/5和2/5频率的信号滤波效果最好。所以根据这个特性,一方面我们要好好利用,一方面也要避免其影响。

信号处理(1) --常用信号平滑去噪的方法

举个应用的例子,比如你的采样频率为10Hz,采样3点滑动平均滤波,但是你的信号在3.3hz左右,你就会发现你的信号被过滤掉了,只留下了噪声。

反之,以美国近期的疫情为例,疫情的采样频率为1天一采样,而且显示出很强的7日一周期的特性(也要过周末信号处理(1) --常用信号平滑去噪的方法)。所以,为了消除这个归一化频率为1/7的噪声影响,采样7点的滑动平均滤波。可以看到所有以7天为一变化的信号分量全部被消除掉了。(下面这个图经常被引用,但是很少有人思考为什么用7天平均方法来平滑数据。)

信号处理(1) --常用信号平滑去噪的方法

回到原本的幅频特性问题上。当点数较少的时候,比如3个点,高频滤波效果并不是很好。所以当取的点数比较少的时候,需要平滑完一遍之后再平滑一遍,直到满意为止。多次平滑之后,高频的衰减非常明显。这也就是说,即使只有3个点平均,多次平滑之后也可以等效为一个较好的低通滤波器。

所以总结一下,移动平均滤波拥有保低频滤高频的特点,而且对于特点频率的滤波具有良好的效果。但是缺点是所有频率分量的信号都会有不同程度衰减。


6、时域和频域的转换关系

时域上的滤波和频域上的滤波是可以互相转换,且一一对应的。也就是时域上的卷积等于频域上的乘积。

下图为3点移动平均滤波法,时域和频域的转换关系:

信号处理(1) --常用信号平滑去噪的方法

虽然前面的 movmean()或者conv()等函数都是用时域实现的信号滤波,但是同样也可以完全在频域上实现。采用ifft(fft(x).*fft(F))实现的滤波效果,和完全时域上的滤波效果是等价的。

信号处理(1) --常用信号平滑去噪的方法信号处理(1) --常用信号平滑去噪的方法信号处理(1) --常用信号平滑去噪的方法

发表评论

相关文章