手把手教系列之梳狀濾波器設(shè)計(jì)實(shí)現(xiàn)
掃描二維碼
隨時(shí)隨地手機(jī)看文章
[導(dǎo)讀]：前面一篇文章關(guān)于IIR/移動(dòng)平均濾波器設(shè)計(jì)的文章。本文來(lái)聊一聊陷波濾波器，該濾波器在混入諧波干擾時(shí)非常有用，算法簡(jiǎn)單，實(shí)現(xiàn)代價(jià)低。本文來(lái)一探其在機(jī)理、應(yīng)用場(chǎng)景。
注：盡量在每篇文章寫(xiě)寫(xiě)摘要，方便閱讀。信息時(shí)代，大家時(shí)間都很寶貴，如此亦可節(jié)約粉絲們的寶貴時(shí)間。
前文所說(shuō)學(xué)習(xí)的倡導(dǎo)2W1H原則，思來(lái)想來(lái)并不全面，本文決定從What Why Where When How幾個(gè)維度展開(kāi)。我稱之為4W1H學(xué)習(xí)法，借鑒管理學(xué)領(lǐng)域中的5W1H，起源于1932年，美國(guó)政治學(xué)家拉斯維爾提出“5W分析法”，后延伸出5W1H法。有興趣的可以找來(lái)閱讀，題外話技術(shù)人員讀一些方法論管理學(xué)方面的書(shū)籍于做人做事個(gè)人認(rèn)為是非常有益的。
梳狀濾波器之What？
在信號(hào)處理中，梳狀濾波器是通過(guò)向其自身添加信號(hào)的延遲而實(shí)現(xiàn)的，從而造成增強(qiáng)或削弱干擾的濾波器。梳狀濾波器的頻率響應(yīng)由一系列規(guī)則間隔的凹口組成，從而呈現(xiàn)出梳狀外觀。其大體拓?fù)湫问饺缦拢?/p>
梳狀濾波器有著大量不同形式的傳遞函數(shù)，其作用是對(duì)周期性信號(hào)增強(qiáng)或削弱周期性信號(hào)，本文主要介紹其中一種形式的Z傳遞函數(shù)
其中：
其信號(hào)流圖如下：
梳狀濾波器英文稱為comb(梳子) filter，這個(gè)名字真是無(wú)與倫比的絕！為何？談到濾波器一定會(huì)重點(diǎn)關(guān)注其對(duì)幅頻響應(yīng)曲線，梳狀濾波器，正是描述其幅頻響應(yīng)的。而幅頻響應(yīng)從本質(zhì)上講是描述系統(tǒng)各頻率能量的放大或者衰減。本文中談到的濾波器就是一個(gè)系統(tǒng)，對(duì)其輸入能量按頻率不同進(jìn)行放大或者衰減，從而起到過(guò)濾作用。
梳狀濾波器之Why？
前面說(shuō)到梳狀濾波器其幅頻響應(yīng)樣子和梳子長(zhǎng)的很像，為啥長(zhǎng)的像，來(lái)一探究竟：
其頻率響應(yīng)為：
現(xiàn)以采樣率20000Hz，10階，阻帶帶寬50Hz為例。其幅頻響應(yīng)曲線如下：
相頻響應(yīng)曲線為：
從幅頻響應(yīng)曲線可看出，其形狀真是如梳子形狀，當(dāng)階數(shù)越大，其齒數(shù)越多。這種形式的梳狀濾波器對(duì)于2KHz，4KHz，6KHz，8KHz，10KHz信號(hào)是衰減作用，也即阻止該頻率通過(guò)，阻帶寬度為50Hz。前面談到了我們可以對(duì)某些頻率信號(hào)加強(qiáng)，而其他的信號(hào)衰減吸收。這里引申出其互補(bǔ)濾波器，其Z傳遞函數(shù)剛好有一點(diǎn)形式上的對(duì)稱性：
其中b為：
此互補(bǔ)濾波器其幅頻響應(yīng)曲線為：
這兩個(gè)濾波器其幅頻響應(yīng)曲線剛好是互補(bǔ)對(duì)稱的。至此，從原理上已基本明了為什么稱其為梳狀濾波器。
梳狀濾波器就其本質(zhì)也是一種IIR濾波器，因?yàn)闉V波器的輸出反饋參與了濾波運(yùn)算。
梳狀濾波器之Where？
費(fèi)這么大勁研究梳狀濾波器，須的知道在什么地方可以去使用它，事實(shí)上梳狀濾波器有著大量的應(yīng)用場(chǎng)合：
-
級(jí)聯(lián)積分梳狀（CIC）濾波器，通常用于插值和抽取操作期間的抗混疊，這些操作會(huì)更改離散時(shí)間系統(tǒng)的采樣率。 -
PAL和NTSC電視解碼器的芯片（也可能是軟件濾波）實(shí)現(xiàn)的2D和3D梳狀濾波器用以降低網(wǎng)點(diǎn)偽影干擾。 -
音頻信號(hào)處理，包括延遲，鑲邊和數(shù)字波導(dǎo)合成中大量應(yīng)用。例如，如果將延遲設(shè)置為幾毫秒，則可以使用梳狀濾波器對(duì)圓柱空腔或振動(dòng)弦中的聲駐波的影響進(jìn)行建模。 -
在天文學(xué)中，天體梳濾波器有望將現(xiàn)有光譜儀的精度提高近一百倍。 -
在聲學(xué)方面，梳齒濾波會(huì)以某些不需要的方式出現(xiàn)。例如，當(dāng)兩個(gè)揚(yáng)聲器在距收聽(tīng)者不同距離處播放同一信號(hào)時(shí)，會(huì)對(duì)信號(hào)產(chǎn)生梳狀濾波效果。在任何封閉的空間中，聽(tīng)眾會(huì)聽(tīng)到直接聲音和反射聲音的混合聲音。由于反射的聲音路徑較長(zhǎng)，因此構(gòu)成了直接聲音的延遲版本，并創(chuàng)建了一個(gè)梳狀濾波器，使兩者在聽(tīng)眾處合并。 -
儀器儀表也常用梳狀濾波器消除諧波干擾，或者選頻放大。 -
......
梳狀濾波器之When？
如遇到諧波干擾，在硬件很難解決時(shí)就可以考慮。
在做圖像處理時(shí)也可以考慮使用該濾波器進(jìn)行精準(zhǔn)重建。
普通單片機(jī)系統(tǒng)
DSP儀器儀表系統(tǒng)中
接下來(lái)就是究竟怎么用C語(yǔ)言以及MATLAB實(shí)現(xiàn)的干貨了，請(qǐng)繼續(xù)......
梳狀濾波器之How？
本文依然借助于matlab的fdatool進(jìn)行設(shè)計(jì)示例：
將其重要設(shè)置標(biāo)注如上，其他的與IIR一文類似，就不羅嗦舉例了。將重要參數(shù)輸入，點(diǎn)擊設(shè)計(jì)就輕松設(shè)計(jì)出相應(yīng)的濾波器參數(shù)了。這里以1000Hz采樣率，40Hz帶寬，20階為例，設(shè)計(jì)出濾波器參數(shù)如下：
% Generated by MATLAB(R) 8.4 and the Signal Processing Toolbox 6.22.
% Generated on: 05-Apr-2020 13:40:29
% Coefficient Format: Decimal
% Discrete-Time IIR Filter (real)
% -------------------------------
% Filter Structure : Direct-Form II
% Numerator Length : 21
% Denominator Length : 21
% Stable : Yes
% Linear Phase : No
Numerator:
0.38970091175151578
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
-0.38970091175151578
Denominator:
1
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0.22059817649696845
C語(yǔ)言實(shí)現(xiàn)非常簡(jiǎn)單，由于梳狀濾波器本質(zhì)上是IIR濾波器，所以也可以直接利用前文ARM的庫(kù)函數(shù)實(shí)現(xiàn)部署。由前面Z傳遞函數(shù)，很容易推導(dǎo)出其差分方程如下：
其互補(bǔ)濾波器差分方程為：
依據(jù)上面公式非常容易設(shè)計(jì)出C代碼，這里將浮點(diǎn)數(shù)實(shí)現(xiàn)共享，如需定點(diǎn)實(shí)現(xiàn)也非常容易，代碼如下：
#include <stdio.h>
#include <math.h>
#include <string.h>
/*長(zhǎng)度應(yīng)為階數(shù)+1*/
#define CMF_RANK 20
typedef double E_SAMPLE;
typedef enum _E_CMF_TYPE{
CMF_TYPE_STOP_BANDS,
CMF_TYPE_PASS_BANDS
}E_CMF_TYPE;
/*定義移動(dòng)平均寄存器歷史狀態(tài)*/
typedef struct _t_CMF
{
E_SAMPLE x[CMF_RANK];
E_SAMPLE y[CMF_RANK];
double b;
double r;
E_CMF_TYPE type;
int index;
}t_CMF;
void comb_filter_init(t_CMF * pCmf,double rho,E_CMF_TYPE type)
{
memset(pCmf,0,sizeof(t_CMF));
pCmf->index = CMF_RANK-1;
pCmf->r = rho;
pCmf->type = type;
if(type==CMF_TYPE_STOP_BANDS)
pCmf->b = (1+rho)/2;
else
pCmf->b = (1-rho)/2;
}
E_SAMPLE comb_filter(t_CMF * pCmf,E_SAMPLE xn)
{
E_SAMPLE yn=0;
int n_N;
int i=0;
n_N = pCmf->index;
if(pCmf->type == CMF_TYPE_STOP_BANDS)
{
/*y[n] = bx[n]-bx[n-N]+ry[n-N]*/
yn = pCmf->b*(xn-pCmf->x[n_N])+pCmf->r*pCmf->y[n_N];
}
else
{
/*y[n] = bx[n]+bx[n-N]+ry[n-N]*/
yn = pCmf->b*(xn+pCmf->x[n_N])+pCmf->r*pCmf->y[n_N];
}
/*存儲(chǔ)yn為下次迭代準(zhǔn)備*/
pCmf->y[n_N] = yn;
pCmf->x[n_N] = xn;
if(pCmf->index==0)
pCmf->index = CMF_RANK-1;
else
pCmf->index--;
return yn;
}
#define SAMPLE_RATE 1000.0f
#define SAMPLE_SIZE 512
#define PI 3.415926f
int main()
{
E_SAMPLE rawSin[SAMPLE_SIZE];
E_SAMPLE outSin[SAMPLE_SIZE];
t_CMF cmf;
FILE *pFile=fopen("./simulationSin.csv","wt+");
if(pFile==NULL)
{
printf("simulationSin.csv opened failed");
return -1;
}
for(int i=0;i<SAMPLE_SIZE;i++)
{
rawSin[i] = 10*sin(2*PI*25*i/SAMPLE_RATE);//+rand()%10;
rawSin[i] += 10*sin(2*PI*50*i/SAMPLE_RATE);
}
/*初始化*/
comb_filter_init(&cmf,-0.22059817649696845,CMF_TYPE_STOP_BANDS);
/*濾波*/
for(int i=0;i<SAMPLE_SIZE;i++)
{
outSin[i]=comb_filter(&cmf,rawSin[i]);
}
for(int i=0;i<SAMPLE_SIZE;i++)
{
fprintf(pFile,"%f,",rawSin[i]);
}
fprintf(pFile,"\n");
for(int i=0;i<SAMPLE_SIZE;i++)
{
fprintf(pFile,"%f,",outSin[i]);
}
/*初始化*/
comb_filter_init(&cmf,-0.22059817649696845,CMF_TYPE_PASS_BANDS);
/*濾波*/
for(int i=0;i<SAMPLE_SIZE;i++)
{
outSin[i]=comb_filter(&cmf,rawSin[i]);
}
/*存儲(chǔ)數(shù)據(jù)*/
fprintf(pFile,"\n");
for(int i=0;i<SAMPLE_SIZE;i++)
{
fprintf(pFile,"%f,",outSin[i]);
}
fclose(pFile);
return 0;
}
同樣利用excel,生成波形如下：
可見(jiàn),經(jīng)過(guò)梳狀濾波器過(guò)濾后，50Hz諧波已經(jīng)被過(guò)濾掉，25Hz 保留下來(lái)，而經(jīng)過(guò)其互補(bǔ)濾波器后，25Hz被過(guò)濾，其50Hz被保留。如看時(shí)域波形不直觀，可利用Python代碼進(jìn)行FFT展開(kāi)分析：
梳妝濾波后FFT譜線圖：
互補(bǔ)梳狀濾波器過(guò)濾后FFT譜線圖：
總結(jié)：
-
梳妝濾波器本質(zhì)上是一種IIR濾波器 -
梳妝濾波器在濾除諧波上，利用極小代價(jià)就可以比較好的濾除諧波干擾 -
其互補(bǔ)濾波器在時(shí)域時(shí)會(huì)失真，使用時(shí)需要考慮 -
如需要考慮計(jì)算效率，可以考慮定點(diǎn)優(yōu)化，但精度會(huì)下降。 -
梳妝濾波器在2D/3D圖像處理廣為應(yīng)用，如有興趣可深入研究 -
如需FFT的python代碼，加作者微信領(lǐng)取
—END—
如果喜歡右下點(diǎn)個(gè)在看，也會(huì)讓我倍感鼓舞
關(guān)注置頂：掃描左下二維碼關(guān)注公眾號(hào)加星
關(guān)注 |
加群 |
免責(zé)聲明：本文內(nèi)容由21ic獲得授權(quán)后發(fā)布，版權(quán)歸原作者所有，本平臺(tái)僅提供信息存儲(chǔ)服務(wù)。文章僅代表作者個(gè)人觀點(diǎn)，不代表本平臺(tái)立場(chǎng)，如有問(wèn)題，請(qǐng)聯(lián)系我們，謝謝！
下載文檔
