手把手教系列之梳狀濾波器設計實現

来源:https://www.cnblogs.com/embInn/archive/2020/05/23/12940708.html

[導讀]:前面一篇文章關於IIR/移動平均濾波器設計的文章。本文來聊一聊陷波濾波器,該濾波器在混入諧波干擾時非常有用,演算法簡單,實現代價低。本文來一探其在機理、應用場景。 註:儘量在每篇文章寫寫摘要,方便閱讀。信息時代,大家時間都很寶貴,如此亦可節約粉絲們的寶貴時間。 前文所說學習的倡導2W1H原則 ...


[導讀]:前面一篇文章關於IIR/移動平均濾波器設計的文章。本文來聊一聊陷波濾波器,該濾波器在混入諧波干擾時非常有用,演算法簡單,實現代價低。本文來一探其在機理、應用場景。

註:儘量在每篇文章寫寫摘要,方便閱讀。信息時代,大家時間都很寶貴,如此亦可節約粉絲們的寶貴時間。

前文所說學習的倡導2W1H原則,思來想來並不全面,本文決定從What Why Where When How幾個維度展開。我稱之為4W1H學習法,借鑒管理學領域中的5W1H,起源於1932年,美國政治學家拉斯維爾提出“5W分析法”,後延伸出5W1H法。有興趣的可以找來閱讀,題外話技術人員讀一些方法論管理學方面的書籍於做人做事個人認為是非常有益的。

梳狀濾波器之What?

在信號處理中,梳狀濾波器是通過向其自身添加信號的延遲而實現的,從而造成增強或削弱干擾的濾波器。 梳狀濾波器的頻率響應由一系列規則間隔的凹口組成,從而呈現出梳狀外觀。其大體拓撲形式如下:

梳狀濾波器有著大量不同形式的傳遞函數,其作用是對周期性信號增強或削弱周期性信號,本文主要介紹其中一種形式的Z傳遞函數\(H(Z)=b\frac{1-Z^{-N}}{1-{\rho}N^{-N}}\)

其中:\(b=\frac{1+\rho}{2}\)

其信號流圖如下:

梳狀濾波器英文稱為comb(梳子) filter,這個名字真是無與倫比的絕!為何?談到濾波器一定會重點關註其對幅頻響應曲線,梳狀濾波器,正是描述其幅頻響應的。而幅頻響應從本質上講是描述系統各頻率能量的放大或者衰減。本文中談到的濾波器就是一個系統,對其輸入能量按頻率不同進行放大或者衰減,從而起到過濾作用。

梳狀濾波器之Why?

前面說到梳狀濾波器其幅頻響應樣子和梳子長的很像,為啥長的像,來一探究竟:

其頻率響應為:

\[H(e^{j\omega})=b\frac{1-e^{j\omega{N}}}{1-\rho{e^{j\omega{N}}}} \]

現以採樣率20000Hz,10階,阻帶帶寬50Hz為例。其幅頻響應曲線如下:

相頻響應曲線為:

![](E:\blog\embInn\DSP\comb Filter\pic\phase12.png)

從幅頻響應曲線可看出,其形狀真是如梳子形狀,當階數越大,其齒數越多。這種形式的梳狀濾波器對於2KHz,4KHz,6KHz,8KHz,10KHz信號是衰減作用,也即阻止該頻率通過,阻帶寬度為50Hz。前面談到了我們可以對某些頻率信號加強,而其他的信號衰減吸收。這裡引申出其互補濾波器,其Z傳遞函數剛好有一點形式上的對稱性:

\[H(Z)=b\frac{1+Z^{-N}}{1-{\rho}N^{-N}} \]

其中b為:

\[b=\frac{1-\rho}{2} \]

此互補濾波器其幅頻響應曲線為:

這兩個濾波器從其幅頻響應曲線剛好是互補對稱的。至此,從原理上已基本明瞭為什麼稱其為梳狀濾波器。

梳狀濾波器就其本質也是一種IIR濾波器,因為濾波器的輸出反饋參與了濾波運算。

梳狀濾波器之Where?

費這麼大勁研究梳狀濾波器,須的知道在什麼地方可以去使用它,事實上梳狀濾波器有著大量的應用場合:

  • 級聯積分梳狀(CIC)濾波器,通常用於插值和抽取操作期間的抗混疊,這些操作會更改離散時間系統的採樣率。
  • PAL和NTSC電視解碼器的晶元(也可能是軟體濾波)實現的2D和3D梳狀濾波器用以降低網點偽影干擾。
  • 音頻信號處理,包括延遲,鑲邊和數字波導合成中大量應用。 例如,如果將延遲設置為幾毫秒,則可以使用梳狀濾波器對圓柱空腔或振動弦中的聲駐波的影響進行建模。
  • 在天文學中,天體梳濾波器有望將現有光譜儀的精度提高近一百倍。
  • 在聲學方面,梳齒濾波會以某些不需要的方式出現。 例如,當兩個揚聲器在距收聽者不同距離處播放同一信號時,會對信號產生梳狀濾波效果。 在任何封閉的空間中,聽眾會聽到直接聲音和反射聲音的混合聲音。 由於反射的聲音路徑較長,因此構成了直接聲音的延遲版本,並創建了一個梳狀濾波器,使兩者在聽眾處合併。
  • 儀器儀錶也常用梳狀濾波器消除諧波干擾,或者選頻放大。
  • ......

梳狀濾波器之How?

本文依然藉助於matlab的fdatool進行設計示例:

將其重要設置標註如上,其他的與IIR一文類似,就不羅嗦舉例了。將重要參數輸入,點擊設計就輕鬆設計出相應的濾波器參數了。這裡以1000Hz採樣率,40Hz帶寬,20階為例,設計出濾波器參數如下:

% 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語言實現非常簡單,由於梳狀濾波器本質上是IIR濾波器,所以也可以直接利用前文ARM的庫函數實現部署。由前面Z傳遞函數,很容易推導出其差分方程如下:

\[\begin{align*} y(n)&=bx(n)-bx(n-N)+\rho{y(n-N)}\\&=\frac{1+\rho}{2}[x(n)-x(n-N)]+\rho{y(n-N)} \end{align*} \]

其互補濾波器差分方程為:

\[\begin{align*} y(n)&=bx(n)-bx(n-N)+\rho{y(n-N)}\\&=\frac{1-\rho}{2}[x(n)+x(n-N)]+\rho{y(n-N)} \end{align*} \]

依據上面公式非常容易設計出C代碼,這裡將浮點數實現共用,如需定點實現也非常容易,代碼如下:

#include <stdio.h>
#include <math.h>
#include <string.h>

/*長度應為階數+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;
/*定義移動平均寄存器歷史狀態*/
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為下次迭代準備*/
    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]);
    }

    /*存儲數據*/
    fprintf(pFile,"\n");
    for(int i=0;i<SAMPLE_SIZE;i++)
    {
        fprintf(pFile,"%f,",outSin[i]);
    }
    fclose(pFile);

    return 0;
}

同樣利用excel,生成波形如下:

可見,經過梳狀濾波器過濾後,50Hz諧波已經被過濾掉,25Hz 保留下來,而經過其互補濾波器後,25Hz被過濾,其50Hz被保留。如看時域波形不直觀,可利用Python代碼進行FFT展開分析:

梳妝濾波後FFT譜線圖:

互補梳狀濾波器過濾後FFT譜線圖:

總結:

  • 梳妝濾波器本質上是一種IIR濾波器
  • 梳妝濾波器在濾除諧波上,利用極小代價就可以比較好的濾除諧波干擾
  • 其互補濾波器在時域時會失真,使用時需要考慮
  • 如需要考慮計算效率,可以考慮定點優化,但精度會下降。
  • 梳妝濾波器在2D/3D圖像處理廣為應用,如有興趣可深入研究

文章出自微信公眾號:嵌入式客棧,更多內容,請關註本人公眾號,嚴禁商業使用,違法必究


您的分享是我們最大的動力!

更多相關文章
  • 為了方便大家寫代碼,vue.js給大家提供了很多方便的修飾符,比如我們經常用到的取消冒泡,阻止預設事件等等~。 目錄 表單修飾符 事件修飾符 滑鼠按鍵修飾符 鍵值修飾符 v-bind修飾符(實在不知道叫啥名字) 表單修飾符 填寫表單,最常用的是什麼?input!v-model~而我們的修飾符正是為了 ...
  • 1、初始化:選中元素進行初始化;$('input').iCheck({ checkboxClass: 'icheckbox_flat-blue', //選擇框的風格。 radioClass: 'icheckbox_flat-blue', });具體顏色:Black — minimal.css //黑 ...
  • 說句心裡話,小編在從事碼農之前。就是乾金融的!而且還是三年“萬金油” 多金多油 金融行業工資高? 我們當初在金融行業做業務,常常會說這樣一句安慰話:三月不出單,一單頂半年。 可想而知,這個薪水是真的讓人眼紅,高薪程式員瞬間就黯淡了。 但我為什麼又轉行做程式員? 一:金融行業就是從去年315開始下滑了 ...
  • 前提: SpringBoot + Vue + ElementUI 實現後臺管理系統模板 -- 前端篇(一):搭建基本環境:https://www.cnblogs.com/l-y-h/p/12930895.html 一、定義公共組件頁面 簡單的頁面效果如下所示:(做的比較粗糙,大致理解頁面即可) 1、 ...
  • 轉載請註明出處:葡萄城官網,葡萄城為開發者提供專業的開發工具、解決方案和服務,賦能開發者。 原文出處:https://blog.bitsrc.io/what-is-deno-and-will-it-replace-nodejs-a13aa1734a74 Deno是什麼? Deno v1.0.0已於5 ...
  • 一、JML初探 ​ 作為一種形式化語言,可以約束 代碼中類和方法的狀態和行為形成規格,通過將一系列具體代碼實現抽象成明確的行為介面,可以形成一種契約式編程模式, 設計者無需考慮實際的數據結構與演算法,可以聚焦於程式的整體邏輯, 形式化語言的無二義性能讓實現者準確理解介面功能,根據問題需要選擇合適的實現 ...
  • 首先安裝Erlang環境 因為 RabbitMQ 需要 erlang 環境的⽀持,所以必須先安裝 erlang 。 如果只是使用RabbitMQ,個人推薦使用RabbitMQ公司維護的 "erlang" 版本,該版本只保留了與RabbltMQ相關的功能, centOS6與7版本的都有,還有erlan ...
  • 【導讀】:前面的文章介紹了移動平均濾波器、IIR濾波器、梳狀濾波器,今天來談談FIR濾波器的設計實現。 本篇文章依然採用4W1H進行描述,從 What Why Where When How 幾個維度展開。為了便於理解4W1H,依然把5W1H的圖附上。 FIR濾波器之What? LTI線性時不變系統沖 ...
一周排行
  • 一:背景 1. 講故事 曾今在項目中發現有同事自定義結構體的時候,居然沒有重寫Equals方法,比如下麵這段代碼: static void Main(string[] args) { var list = Enumerable.Range(0, 1000).Select(m => new Point ...
  • 最近一個朋友有個關於素數的小東西要寫一下,素數是什麼呢?除了1和他本身不能被其他數整除,那麼這個數就是素數,1除外哦。我們知道概念那就很簡單了,直接代碼擼起。 ...
  • 前言 在開發編程中,我們經常會遇到功能非常相似的功能模塊,只是他們的處理的數據不一樣,所以我們會分別採用多個方法來處理不同的數據類型。但是這個時候,我們就會想一個問題,有沒有辦法實現利用同一個方法來傳遞不同種類型的參數呢? 這個時候,泛型也就因運而生,專門來解決這個問題的。 泛型是在C 2.0就推出 ...
  • 本文章主要用於介紹在Asp.Net Mvc(C#)中使用Fleck製作一個Html5的即時聊天室,含有完整代碼和演示Demo。 ...
  • 出庫單的功能。能學習了出庫單管理之後,WMS的 主體功能算是完成了。當然一個成熟的WMS還包括了盤點,報表,策略規則,移庫功能及與其他系統(ERP、TMS等)的介面,實現無縫集成,打破信息孤島,讓數據實時、準確和同步。 ...
  • Data StructureThere're two types of variables in C#, reference type and value type.Enum:enum Color{Red=0,Green=1}//equals to enum Color{Red,//start fr... ...
  • 0. 前言 該項目使用Maven進行管理和構建,所以需要預先配置好Maven。嗯,在這個系列里就不做過多的介紹了。 1. 創建項目 先創建一個pom.xml 文件,添加以下內容: <?xml version="1.0" encoding="UTF-8"?> <project xmlns="http: ...
  • API 概述 API(Application Programming Interface),應用程式編程介面。 Java API是一本程式員的 字典 ,是JDK中提供給我們使用的類的說明文檔。 這些類將底層的代碼實現封裝了起來,我們不需要關心這些類是如何實現的,只需要學習這些類如何使用即可。 所以我 ...
  • 女程式員是這麼徵婚的: SELECT * FROM 男人們 WHERE 未婚=true and 同性戀=false and 有房=true and 有車=true and 條件 in (帥氣,紳士,大度,氣質,智慧,溫柔,體貼,會浪漫,活潑,可愛,最好還能帶孩子) and 年齡 between(24 ...
  • 有很多剛學習軟體測試的小伙伴,都會在網路上找尋各種學習資料,去提升自己的專業技能水平。因此,我決定定期分享我整理收集的一些軟體測試的測試工具下載、面試寶典、視頻教學合集。都整理好了,有需要的可以關註我(獲取方式在文末) 軟體測試的學習,不止是基礎理論,還需要學習測試工具的用法,如介面工具Postma ...