fir1截止频率计算_如何快速设计一个FIR滤波器(二)
一、理想低通濾波器單位脈沖響應是什么樣
在如何快速設計一個FIR濾波器(一)中,我們介紹了一種簡單設計FIR的方法——零極點法。這個方法非常簡單,稍加培訓,用筆和紙就能完成;當然缺點也很顯而易見:零極點設計出的濾波器,只能給出大概的頻率響應,對于一些要求較高的系統,顯得無能為力。今天我們介紹一種更加嚴謹的方法。
現在假設我們要設計一個低通濾波器,截止頻率為
,其理想頻率響應可以用如下函數表示:
則該濾波器的脈沖響應為:
可見脈沖響應為一個sinc函數。畫出來大概張這個樣子:
這個函數非常重要哦,建議都自己手動推導一下,非常有意思,我們把幅值最大的那一瓣稱之為主瓣,剩下稱之為從瓣,是不是很形象呢?注意主瓣的周期是從瓣的2倍哦。
可能你已經看見了,脈沖響應是由無窮多個——不對呀,我們是想設計一個有限脈沖響應的,結果出現了無窮多個響應,這可如何搞?——既然臣妾做不到把所有的脈沖響應都用上(因為我是有限脈沖濾波器FIR啊),那我只能截取有限一部分進行代替了,也就是——加窗。
窗函數這個東西很不好理解,我們就多花點功夫說說這個事。
二、什么是加窗
現在假設有一個函數,表達式為
,為簡單起見,令
,
,很顯然,余弦函數的頻域包含兩個分量
,即
Hz。
假如我們現在以10Hz的頻率進行采樣,顯然是滿足香濃采樣定理的,采樣脈沖及其頻譜如下圖:
采樣過程就是拿采樣信號和原始信號進行乘積,那采樣之后的信號長啥樣呢?
上圖即為采樣后的信號及其頻譜。注意:此時采樣后的信號是在時域無限擴展的,頻譜也一樣。但實際的計算中是不可能處理無窮多數量的信號的,那怎么辦呢?——截短唄,只觀察一部分樣本行不?貌似我們也沒有更好的選擇了。
那怎么截短呢?最簡單的方式就是矩形函數了,如下圖所示:
矩形函數的方法簡單粗暴,只保留一部分(幅值為1),其余全部設為0,然后那矩形函數和采樣信號相乘,就得到加窗(矩形窗)之后的采樣信號了。矩形信號的頻譜就是我們熟悉的sinc函數。加窗之后的采樣信號及其頻譜如下:
加窗后的離散余弦信號可以看成窗函數(矩形,頻譜為sinc函數)與余弦離散信號的乘積。那在頻域呢就是sinc函數和余弦頻譜的卷積。看一個圖就明白加窗后的頻譜是怎么來的了:
可以看出,由于加窗行為,頻譜中頻率成分變復雜了,由于矩形窗函數的頻譜是無限寬的,在頻譜的最大值處也出現了一定的畸變(略大于1)。這種有限長的離散信號在時域中計算機是可以處理(離散的),但是在頻域還是不能,因為還是連續信號。在如何理解離散傅里葉變換及Z變換一文中,我們介紹了將有限長的信號進行周期延拓,就得到一個周期的離散信號,其傅里葉變換也是周期離散的,這樣就可以在計算機中處理了。關于離散傅里葉變換,可以參照:J Pan:如何理解離散傅里葉變換及Z變換?zhuanlan.zhihu.com
現在問題就變成了如何進行周期延拓——卷積,什么意思呢?我們現在找一個采樣脈沖,時域中周期和加窗后的采樣信號的周期一致(都是1),其形貌如下圖(左)所示:
我們已經知道,周期延拓采樣脈沖的頻譜仍未采樣脈沖如上圖(右)所示。現在就可以做一個有意思的事情:周期延拓在數學上可以看成是一個有限長的信號(長度即為周期)和一個同樣周期的無限長采樣脈沖的卷積。什么意思呢?看一個圖就一目了然了:
感興趣的可以用卷積的數學定義推導一下看看是不是這樣,很簡單的。由卷積定理我們知道,卷積和乘積在時域和頻域是對偶的,也就是說剛才我們所做的周期延拓(時域卷積)在頻域是乘積哦,具體如下圖所示。
所以最終得到的經過周期延拓的加窗信號及其頻譜如下圖所示。可以看出,進過周期延拓,加窗導致的頻譜泄漏似乎得到了抑制,加窗對頻譜沒什么影響。
事實是不是這樣呢?——直覺告訴我們應該不是。我們之所以最終得到這個結論,是因為我們選的初始函數太簡單了。余弦函數是周期信號,頻譜是有限寬度的,如果窗函數的寬度是信號周期的整數倍,則頻域采樣(時域周期延拓)將正好采在sinc函數的從瓣取值為零的地方,泄漏的頻譜沒有采到而已,具體如下圖所示。
可能心細有童鞋會有疑問,為啥主瓣的最大值會有一定的偏移?——原因就是sinc函數主瓣和從瓣的周期不一樣,卷積時主瓣與從半疊加,最大值就偏移了。
實際工程中,我們可能不知道被處理信號的周期,或者壓根就不是周期的信號,那會出現什么呢?現在假如我們把窗函數加寬,比如增加到原來的1.4倍,如下圖所示:
則可以得到加窗后的信號及頻譜為:
加窗并進行周期延拓后的信號及頻譜為:
可見,由于窗函數寬度和原信號周期不一致,加窗后的信號進行周期延拓時出現了不連續,頻域采樣時不像前面那樣只采集到sinc為零的地方,頻譜出現了較大的泄漏,如下圖所示。
實際工程上絕大多數信號都不是周期信號,也不是有限長信號,而是我們只能觀察有限長度而已。要想不出現頻譜泄漏,需要窗函數長度是被處理信號頻譜所有分量的整數倍,這幾乎是不能達到的。因此加窗后,頻譜泄漏是不可避免的,只能盡量減小。本節中所用的算例來自離散傅里葉變換DFT基本原理圖解,感興趣請前往閱讀。
三、窗函數都有哪些
最開始我們介紹了一個理想的低通濾波器,在頻域內其幅值響應是一個矩形。這個矩形進行傅里葉逆變換,發現單位沖擊響應有無窮多個,我們沒辦法,用了一個矩形窗進行截短,注意著兩個矩形不是同一個事情哦,一個是幅值響應(頻域),一個截短函數(時域),只不過樣子都呈現矩形。
現在我們都放在頻域來看:首先時域矩形函數在頻域為sinc函數,時域的矩形函數截短(乘積)在頻域對應卷積。也就是說在頻域里,脈沖響應截短后的濾波器頻域幅值響應是矩形信號和sinc信號的卷積,仔細想想這段話哦!下圖展示了卷積之后的信號的樣子。
可以想象:假設窗函數寬度無窮大,則對應的sinc函數無窮窄,則卷積之后的函數非常接近理想矩形;反之,若窗函數很窄,在sinc函數就會變得很寬,則卷積之后的函數與理想矩形差別就比較大,這就是傳說中的不確定原理哦,具體可參照:J Pan:如何理解不確定性原理—不確定or測不準??zhuanlan.zhihu.com
那什么樣的窗函數比較好呢?——在頻域看就是要能量相對集中,也就是旁瓣要低。因為出現泄漏的主要原因就是窗函數的頻譜是無限長的,與信號的頻譜卷積時,主瓣與從瓣會出現疊加,因此從瓣的能量越小,影響就越小。在時域看就是加窗后的函數在進行周期延拓時盡量不要有非連續,因為非連續就需要很多分量來模擬,就造成了頻譜泄露。那什么樣函數的頻譜主瓣能量較集中呢?——主要有以下是幾種常見的窗函數:矩形窗、漢寧窗、漢明窗、凱撒窗等。
這些窗函數對應的頻譜分別為:
現在我們拿漢寧窗(Hanning)來舉個例子。還是原來的余弦函數
,
,
。假設窗函數的寬度
,Hanning窗窗函數及其頻譜為:
可見其頻譜中主瓣能量與矩形窗比,集中了很多。余弦信號采樣、加窗之后的信號及其頻譜為:
由于Hanning窗從零開始,結束的時候也是零,因此加窗之后的信號開始和結束也都是零,連續性好了,從頻譜看,泄漏也少了。這個從周期延拓后(DTFT)的信號與頻譜看的更清楚。
四、如何基于窗函數法設計FIR
這篇文章中,我們主要說怎么較為精確的設計一個FIR,我們花了大量篇幅說了窗函數,這是因為窗函數在實際的工程應用中很多,比如基于窗函數的FIR就是FIR設計中最廣泛運用的一個方法。接下來,我們要介紹一下什么是用窗函數設計FIR。
在文章開始的時候,我們已經介紹了,對于理想低通濾波器,其脈沖響應為sinc函數,且有無窮多個,我們必須加窗才能處理,加窗就會有一定的頻譜泄漏,因此實際設計出來的濾波器和理想濾波器還是有差別的。
基于窗函數的FIR設計就是通過選擇合適常函數來找到滿足要求的濾波器,一般步驟是這樣的:
1)確定頻域的響應函數
,低通、高通、帶通或者其他;
2)確定頻域的響應函數
的傅里葉逆變換,找到連續脈沖響應函數
;
3)對脈沖響應函數
按照一定的采樣頻率進行采樣,獲得離散信號
;
4)選擇合適的窗函數,對離散信號
加窗,獲得有限長度的脈沖響應
;
當然實際實踐中,我們只需要知道這個過程,具體細節不需要我們考慮,因為MATLAB提供了一個強大的函數,幫我們都做好了——fir1。
fir1的基本語法如下:h = fir1(n,Wn,ftype,window)
其中n表示濾波器的階數(order),Wn表示歸一化后的濾波器的截止頻率,可表示成
的形式。舉個例子,比如一個帶狀濾波器,采樣頻率
是200Hz,兩個截止頻率分別為10Hz和50Hz,則
,即截止頻率除以
。ftype表示濾波器類型,可以是lowpass, highpass, bandpass, bandstop, or multiband filter。window表示窗函數類型,前面我們列到的窗函數都可以選擇。h為所設計的濾波器的系數。
舉個簡單的例子:
h = fir1(48,[0.35 0.65]);
freqz(h,1,512)
這是一個48階的帶通濾波器,歸一化的截止頻率為[0.35 0.65],其幅值響應和相位響應如下圖所示。
可見,采用fir1函數設計FIR濾波器非常方便。
除了fir1函數,fir2函數也有時候會用到,fir2函數是什么意思呢?——我們稱之為基于頻率采樣法的設計方法,具體做法就是:把濾波器期望的幅值響應用一組分立的點
來表示,
表示頻率為
時的期望幅值響應,然后在不同的點之間進行線性插值,得到完成的期望幅值響應,然后進行傅里葉逆變換并用Hamming(漢明窗)截短來獲得濾波器的系數。
基本語法為:h = fir2(n,f,m)
n表示濾波器階數,f表示分離點的矢量,m表示分立點對應的幅值響應矢量,舉個例子就明白了。現在要設計一個低通濾波器,幅值響應如下圖所示:
顯然四個特征點的坐標分別為(0,1), (0.2,1), (0.2,0), (1,0) ,于是我們可以獲得分立點的矢量為f=[0 0.2 0.2 1] ,對應幅值響應矢量為m=[1 1 0 0]。
f=[0 0.2 0.2 1] ;
m=[1 1 0 0];
h = fir2(48,f,m);
freqz(h,1)
五、如何基于響應最優法設計FIR
前面我們說了兩種基于窗函數的方法,主要是用窗函數對無限長的脈沖進行截短,獲得近似的頻域響應。除了窗函數法,還有另外的一種解決方案也比較常用,我們稱之為“最優法”,主要思路就是找到一組脈沖響應,讓它的頻域響應
與期望的濾波器的頻域響應
盡可能的一致,主要通過兩種方法來實現,一個是最小二乘法,另一個是切比雪夫法。
(一)最小二乘法
與頻率采樣法近似,期望的頻率響應用一組分立的點
來表示,在不同的點之間進行插值,然后來尋找濾波器的系數
時期滿足
其中
為不同頻率下的權重。
MATLAB指令為:h=firls(n,f,m)
n為濾波器階數,f表示分離點的矢量,m表示分立點對應的幅值響應矢量。仍以前面說的低通濾波器為例:
f=[0 0.2 0.2 1] ;
m=[1 1 0 0];
h = firls(48,f,m);
freqz(h,1)
(二)切比雪夫法(又叫帕克斯-麥克勞林法)
與最小二乘法不同(方差最小),切比雪夫法采用的方案是最大誤差最小,即:
MATLAB指令為:h=firpm(n,f,m)
n、f、m定義和前面一致,不同的是頻率響應在同一頻率下必須去不同的值。前面的例子中當頻率為0.2(歸一化后)時,頻率響應可直接從1降到0,對于切比雪夫法,規定不能這樣做,我們將第二個0.2改為0,21。
f=[0 0.2 0.21 1] ;
m=[1 1 0 0];
h = firpm(48,f,m);
freqz(h,1)
六、如何利用MATLAB設計FIR濾波器
要想快速設計一個FIR濾波器,MATLAB可以說是一個最有利的幫手,當你知道較多濾波器知識和相關函數時,可以直接調用函數進行設計,靈活且方便。當你不知道時,也沒有關系,打開MATLAB,在命令行輸入:filterDesigner或者fdatool(老版MATLAB),你就能看到如下界面:
通過勾勾選選就能設計濾波器哦!
如果你又進步了一點,已經設計出一個濾波器,相看看性能怎么樣,也簡單,打開MATLAB,輸入:fvtool(h),h就是你設計的濾波器的系數,你就會看到各種你想要的濾波器特性:幅值響應、相位響應、脈沖響應、零極點等等。
總結
以上是生活随笔為你收集整理的fir1截止频率计算_如何快速设计一个FIR滤波器(二)的全部內容,希望文章能夠幫你解決所遇到的問題。
- 上一篇: Oracle 11g 单机安装 —— 文
- 下一篇: vue中设置html的fontsize,