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

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

关键字：信号；去噪；Matlab

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

1、移动平均法

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

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

以3点平均（窗口长度为3）公式为例，原数据为x，平滑后的数据为y：

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

对y(n)和y(n+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)

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)

中间部分两者完全一致，但是两端有所差别。主要是因为，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 = filter(b,a,x)

其中b和a的定义为：

其中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。

5、移动平均的幅频响应

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

以中心窗口为例，

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)
end
H_iw_Sum = abs(sum(H_iw,1));%相加
figure()
plot(w/2/pi,H_iw_Sum)

由于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频率的信号滤波效果最好。所以根据这个特性，一方面我们要好好利用，一方面也要避免其影响。

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

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

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

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

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

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

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

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

这也意味着你也可以在频域上操作，实现想要的滤波。比如想要低频通过高频衰减，就把fft后的信号，高频部分强行等于0即可。比如想要消除某个频率的信号（陷波），就令fft后那个信号的频率等于0即可。同理，想要把振幅衰减1/2，就在对应频域上乘以0.5.