滑动平均法

姓名：武志霞.      学号：20021110081

转载自：https://blog.csdn.net/weixin_42943114/article/details/107693068，有删节。

【嵌牛导读】滑动平均法是一种信号平滑去噪方法。

【嵌牛鼻子】滑动平均法

【嵌牛提问】滑动平均法原理及应用

【嵌牛正文】：

参考书目：

1《数字信号处理》-胡广书

2《Digital Signal Processing - A Practical Guide For Engineers and Scientists》 - Steven W.Smith

1.0 移动平均法的方法原理

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

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

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

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

当然这两者都是等价的。

1.1 matlab内自带函数实现移动平均法

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

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

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

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

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

y = movmean( x , 5 );

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

%移动平均滤波

clear

clc

close all

N_window = 5;%窗口长度(最好为奇数）

t = 0:0.1:10;

A = cos(2*pi*0.5*t)+0.3*rand(size(t));

B1 = movmean(A,N_window);

figure(1)

plot(t,A,t,B1)

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

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

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

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

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

可以验证：

N_window = 5;%窗口长度(最好为奇数）

t = 0:0.25:10;%时间

A = cos(2*pi*0.5*t)+0.3*rand(size(t));

B1 = movmean(A,N_window);

F_average = 1/N_window*ones(1,N_window);%构建卷积核

B2 = conv(A,F_average,'same');%利用卷积的方法计算

figure(2)

plot(t,B1,t,B2)

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

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

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

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

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/5,1/5,1/5,1/5,1/5] , [1] , x );

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

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

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

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

1.4 移动平均的幅频响应

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

y(n)=1/3（x(n−1)+x(n)+x(n+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)

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个点，高频滤波效果并不是很好。所以当取的点数比较少的时候，需要平滑完一遍之后再平滑一遍，直到满意为止。这个原理也可以通过幅频特性来解释，多次平滑相当于乘了多个H(z)。对于平滑2次，它的图像也就是abs(sum(H_iw,1).*sum(H_iw,1))；对于平滑4次，它的图像相当于乘以四个sum(H_iw,1)。（注：因为时域上的卷积等于频域上的乘积，多次卷积对应着多次乘积。）

可以看到，多次平滑之后，高频的衰减非常明显。这也就是说，即使只有3个点平均，多次平滑之后也可以等效为一个较好的低通滤波器。

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

1.5 时域和频域的转换关系

额外补充一部分小内容，可能前面有些概念加入的太突然。很多人可能觉得之前时域上的平均法非常好理解，为什么突然加入幅频特性图，又是Z变换又是fft的。

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

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

同样的滤波操作，可以用时域公式：y（n）=1/3 (x(n−1)+x(n)+x(n+1))，进行描述。也可以用频域上，滤波器的幅频特性进行描述。

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

下面是展示了窗口长度为3的平滑滤波，从时域上和频域上对信号进行滤波的对比：

%实验，检验频域和时域的一致性

%以3点滤波为例

clear

clc

close all

N_window = 3;%窗口长度(最好为奇数）

t = 0:0.1:10;

A = cos(2*pi*0.3*t)+0.1*cos(2*pi*5*t)+0.2*randn(size(t));

F_average = 1/N_window*ones(1,N_window);%创建滤波器

B2 = conv(A,F_average,'same');%利用卷积的方法计算

figure(1)

plot(t,A,'k','linewidth',0.8)

%计算原信号的fft

A_fft=fft(A);

%构建频域上的滤波器

F_average2=zeros(size(t));%长度与x相同，为了后面.*运算

F_average2(1:(N_window-1)/2+1) = 1/N_window;

F_average2(end-(N_window-1)/2+1:end) = 1/N_window;%前后设置对称，使得相位不变

F_Fft = fft(F_average2);

figure(2)

subplot(1,2,1)

plot(linspace(0,1,length(t)),abs(A_fft),'linewidth',1);

subplot(1,2,2)

plot(linspace(0,1,length(t)),abs(F_Fft),'linewidth',1);

%进行反逆变换

B3=ifft(A_fft.*F_Fft);

figure()

plot( t,B2,t,B3 )%对比时域操作和频域操作的差异

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