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