姓名：武志霞.? ? ? 學(xué)號(hào)：20021110081
轉(zhuǎn)載自：https://blog.csdn.net/weixin_42943114/article/details/107693068，有刪節(jié)。
【嵌牛導(dǎo)讀】滑動(dòng)平均法是一種信號(hào)平滑去噪方法。
【嵌牛鼻子】滑動(dòng)平均法
【嵌牛提問(wèn)】滑動(dòng)平均法原理及應(yīng)用
【嵌牛正文】：
參考書目：
1《數(shù)字信號(hào)處理》-胡廣書
2《Digital Signal Processing - A Practical Guide For Engineers and Scientists》 - Steven W.Smith
1.0 移動(dòng)平均法的方法原理
滑動(dòng)平均法（moving average）也叫做移動(dòng)平均法、平均法、移動(dòng)平均值濾波法等等，是一種時(shí)間域思想上的信號(hào)光滑方法。算法思路為，將該點(diǎn)附近的采樣點(diǎn)做算數(shù)平均，作為這個(gè)點(diǎn)光滑后的值。
一般窗口為對(duì)稱窗口，防止出現(xiàn)相位偏差。窗口一般為奇數(shù)。
以3點(diǎn)平均（窗口長(zhǎng)度為3）公式為例，原數(shù)據(jù)為x，平滑后的數(shù)據(jù)為y：
對(duì)y(n)和y(n+1)相減，可以得到另一種計(jì)算形式：
當(dāng)然這兩者都是等價(jià)的。
1.1 matlab內(nèi)自帶函數(shù)實(shí)現(xiàn)移動(dòng)平均法
matlab有兩個(gè)函數(shù)實(shí)現(xiàn)滑動(dòng)平均法，一個(gè)是smoothdata()函數(shù)，一個(gè)是movmean()函數(shù)。
以窗口長(zhǎng)度為5為例，smoothdata()函數(shù)調(diào)用方法為：
y = smoothdata( x , 'movmean' , 5 );
但是這個(gè)smoothdata函數(shù)實(shí)際上是調(diào)用了movmean()函數(shù)。所以如果直接使用的話，直接用movmean()會(huì)更快。
movmean()函數(shù)的調(diào)用方法為：
y = movmean( x , 5 );
下面以一個(gè)加噪聲的正弦信號(hào)為例：
%移動(dòng)平均濾波
clear
clc
close all
N_window = 5;%窗口長(zhǎng)度(最好為奇數(shù)）
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 利用卷積函數(shù)conv()實(shí)現(xiàn)移動(dòng)平均法
根據(jù)之前的公式，我們可以看到
y(n)=1/3?(x(n?1)+x(n)+x(n+1))
y(n)=1/3?(x(n?1)+x(n)+x(n+1))
就相當(dāng)于一個(gè)對(duì)x和向量[1/3 1/3 1/3]做卷積。
可以驗(yàn)證：
N_window = 5;%窗口長(zhǎng)度(最好為奇數(shù)）
t = 0:0.25:10;%時(shí)間
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);%構(gòu)建卷積核
B2 = conv(A,F_average,'same');%利用卷積的方法計(jì)算
figure(2)
plot(t,B1,t,B2)
中間部分兩者完全一致，但是兩端有所差別。主要是因?yàn)椋琺ovmean()函數(shù)在處理邊緣時(shí)，采用減小窗口的方式，而conv()相當(dāng)于在兩端補(bǔ)零。所以如何處理邊緣也是值得注意的。
1.3 利用filter濾波函數(shù)實(shí)現(xiàn)移動(dòng)平均法
首先介紹一下Z變換。以向前的滑動(dòng)平均為例（這里中間值不是n而是n+1，所以相位會(huì)移動(dòng)）。
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變換可以簡(jiǎn)單的理解為，把x(n+k)替換為z^(-k)，即
因此對(duì)于filter濾波函數(shù)，輸入的格式為：
y = filter(b,a,x)
其中b和a的定義為：
其中a1必須為1。所以對(duì)應(yīng)的移動(dòng)平均濾波可以表示為：
y = filter( [1/5,1/5,1/5,1/5,1/5] , [1] , x );
它和下面代碼的是等價(jià)的（在邊緣上的處理方式有所不同
y = movmean( x , [4,0] );
代表有偏移的滑動(dòng)平均濾波，4是向后4個(gè)點(diǎn)的意思，0是向前0個(gè)點(diǎn)的意思。
因?yàn)?filter濾波器使用有偏移的向后濾波。濾波后，相位會(huì)發(fā)生改變。所以通常采用零相位濾波器進(jìn)行濾波，matlab內(nèi)的函數(shù)為filtfilt()。原理從函數(shù)名字上就可以體現(xiàn)出來(lái)，就是先正常濾波，之后再將信號(hào)從后向前再次濾波。正濾一遍反濾一遍，使得相位偏移等于0。
1.4 移動(dòng)平均的幅頻響應(yīng)
幅頻響應(yīng)可以通過(guò)之前1.3得到的H(z)函數(shù)來(lái)得到，在單位圓上采樣，也就是把z替換為e^iw。以中心窗口為例，
y(n)=1/3（x(n?1)+x(n)+x(n+1))
H(iw)的絕對(duì)值就是該濾波方法的幅頻響應(yīng)。以3點(diǎn)濾波為例，matlab代碼為：
%計(jì)算頻率響應(yīng)
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變換在單位圓上的特性相當(dāng)于傅里葉變換，所以上面代碼等效于下面計(jì)算方法：
%計(jì)算頻率響應(yīng)
N_window = 3;
Y=zeros(400,1);
Y(1:N_window) = 1/3;%設(shè)置濾波器
Y_F = fft(Y);
figure()
plot(linspace(0,0.5,200),abs(Y_F(1:200)));
matlab也有自帶的函數(shù)來(lái)看頻率特性，freqz()，推薦使用這種。
其中，歸一化頻率等于信號(hào)頻率除以采樣頻率f/Fs，采樣頻率等于時(shí)間采樣間隔的倒數(shù)1/dt。對(duì)比不同窗口長(zhǎng)度的幅頻響應(yīng)，可以看到：
1平均所采用的點(diǎn)數(shù)越多，高頻信號(hào)的濾波效果越好。
2 3點(diǎn)平均對(duì)于1/3頻率的信號(hào)濾波效果最好，5點(diǎn)平均對(duì)1/5和2/5頻率的信號(hào)濾波效果最好。所以根據(jù)這個(gè)特性，一方面我們要好好利用，一方面也要避免其影響。
舉個(gè)應(yīng)用的例子，比如你的采樣頻率為10hz，采樣3點(diǎn)滑動(dòng)平均濾波，但是你的信號(hào)在3.3hz左右，你就會(huì)發(fā)現(xiàn)你的信號(hào)被過(guò)濾掉了，只留下了噪聲。
反之，以美國(guó)近期的疫情為例，疫情的采樣頻率為1天一采樣，而且顯示出很強(qiáng)的7日一周期的特性（病毒也要過(guò)周末）。所以，為了消除這個(gè)歸一化頻率為1/7的噪聲影響，采樣7點(diǎn)的滑動(dòng)平均濾波。可以看到所有以7天為一變化的信號(hào)分量全部被消除掉了。
（下面這個(gè)圖經(jīng)常被引用，但是很少有人思考為什么用7天平均方法來(lái)平滑數(shù)據(jù)。）
回到原本的幅頻特性問(wèn)題上。當(dāng)點(diǎn)數(shù)較少的時(shí)候，比如3個(gè)點(diǎn)，高頻濾波效果并不是很好。所以當(dāng)取的點(diǎn)數(shù)比較少的時(shí)候，需要平滑完一遍之后再平滑一遍，直到滿意為止。這個(gè)原理也可以通過(guò)幅頻特性來(lái)解釋，多次平滑相當(dāng)于乘了多個(gè)H(z)。對(duì)于平滑2次，它的圖像也就是abs(sum(H_iw,1).*sum(H_iw,1))；對(duì)于平滑4次，它的圖像相當(dāng)于乘以四個(gè)sum(H_iw,1)。（注：因?yàn)闀r(shí)域上的卷積等于頻域上的乘積，多次卷積對(duì)應(yīng)著多次乘積。）
可以看到，多次平滑之后，高頻的衰減非常明顯。這也就是說(shuō)，即使只有3個(gè)點(diǎn)平均，多次平滑之后也可以等效為一個(gè)較好的低通濾波器。
所以總結(jié)一下，移動(dòng)平均濾波擁有保低頻濾高頻的特點(diǎn)，而且對(duì)于特點(diǎn)頻率的濾波具有良好的效果。但是缺點(diǎn)是所有頻率分量的信號(hào)都會(huì)有不同程度衰減。
1.5 時(shí)域和頻域的轉(zhuǎn)換關(guān)系
額外補(bǔ)充一部分小內(nèi)容，可能前面有些概念加入的太突然。很多人可能覺得之前時(shí)域上的平均法非常好理解，為什么突然加入幅頻特性圖，又是Z變換又是fft的。
其實(shí)時(shí)域上的濾波和頻域上的濾波是可以互相轉(zhuǎn)換，且一一對(duì)應(yīng)的。也就是時(shí)域上的卷積等于頻域上的乘積。
下圖為3點(diǎn)移動(dòng)平均濾波法，時(shí)域和頻域的轉(zhuǎn)換關(guān)系：
同樣的濾波操作，可以用時(shí)域公式：y（n）=1/3 (x(n?1)+x(n)+x(n+1))，進(jìn)行描述。也可以用頻域上，濾波器的幅頻特性進(jìn)行描述。
雖然前面的 movmean()或者conv()等函數(shù)都是用時(shí)域?qū)崿F(xiàn)的信號(hào)濾波，但是同樣也可以完全在頻域上實(shí)現(xiàn)。采用ifft(fft(x).*fft(F))實(shí)現(xiàn)的濾波效果，和完全時(shí)域上的濾波效果是等價(jià)的。
下面是展示了窗口長(zhǎng)度為3的平滑濾波，從時(shí)域上和頻域上對(duì)信號(hào)進(jìn)行濾波的對(duì)比：
%實(shí)驗(yàn)，檢驗(yàn)頻域和時(shí)域的一致性
%以3點(diǎn)濾波為例
clear
clc
close all
N_window = 3;%窗口長(zhǎng)度(最好為奇數(shù)）
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);%創(chuàng)建濾波器
B2 = conv(A,F_average,'same');%利用卷積的方法計(jì)算
figure(1)
plot(t,A,'k','linewidth',0.8)
%計(jì)算原信號(hào)的fft
A_fft=fft(A);
%構(gòu)建頻域上的濾波器
F_average2=zeros(size(t));%長(zhǎng)度與x相同，為了后面.*運(yùn)算
F_average2(1:(N_window-1)/2+1) = 1/N_window;
F_average2(end-(N_window-1)/2+1:end) = 1/N_window;%前后設(shè)置對(duì)稱，使得相位不變
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);
%進(jìn)行反逆變換
B3=ifft(A_fft.*F_Fft);
figure()
plot( t,B2,t,B3 )%對(duì)比時(shí)域操作和頻域操作的差異
這也意味著你也可以在頻域上操作，實(shí)現(xiàn)想要的濾波。比如想要低頻通過(guò)高頻衰減，就把fft后的信號(hào)，高頻部分強(qiáng)行等于0即可。比如想要消除某個(gè)頻率的信號(hào)（陷波），就令fft后那個(gè)信號(hào)的頻率等于0即可。同理，想要把振幅衰減1/2，就在對(duì)應(yīng)頻域上乘以0.5。
