滑动平均法
姓名：武志霞.学号：***********
【嵌⽜导读】滑动平均法是⼀种信号平滑去噪⽅法。
【嵌⽜⿐⼦】滑动平均法
【嵌⽜提问】滑动平均法原理及应⽤
【嵌⽜正⽂】：
参考书⽬：
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。
linspace函数调用的格式为
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天平均⽅法来平滑数据。）

版权声明：本站内容均来自互联网，仅供演示用，请勿用于商业和其他非法用途。如果侵犯了您的权益请与我们联系QQ:729038198，我们将在24小时内删除。