FIR濾波器，低通、高通、帶通、帶阻VC實現(轉)

1.前言：數字信號處理相關知識准備

通常來說，一種理想濾波器的頻率響應是很容易理解的，如圖所示。
 

 

 

圖1 濾波器頻響
以低通為例，濾波器頻率響應函數為

所謂濾波器處理的過程，簡單來說，可以用公式

 

來表示，由卷積的性質可以知道，該公式的另一種形式為

其中x(n)為要處理的數據序列，h(n)為逼近濾波器的時域響應

其中，hd(n)為對應不同類型濾波器的單位沖擊響應，比如說低通的hd(n)為sinc函數。
我們知道，高通可以有全通減低通得到，帶通可由兩個低通相減得到，帶阻可由低通加高通得到。
————————————————

2. 具體VC實現過程
有了上面簡單的回顧之后，我們就可以進行VC上濾波器的實現了。首先是hd(n)的實現，具體代碼如下:
頭文件聲明部分

#pragma once
class CFIRWIN
{
public:
 CFIRWIN(void);
 ~CFIRWIN(void);
 void firwin(int n,int band,int wn,double h[],double kaiser=0.0,double fln=0.0,double fhn=0.0);
 double window(int type,int n,int i,double beta);//窗函數的計算
 double kaiser(int i,int n,double beta);
 double bessel0(double x);
};

 

源文件實現部分

#include "StdAfx.h"

#include "FIRWIN.h"

#include <math.h>

 

CFIRWIN::CFIRWIN(void)

{

}

 

 

CFIRWIN::~CFIRWIN(void)

{

}

 

void CFIRWIN::firwin(int n,int band,int wn,double h[],double kaiser,double fln,double fhn)

{

    int i,n2,mid;

    double s,pi,wc1,wc2,beta,delay,fs;

    fs=44100;//44kHz

    beta=kaiser;

    pi=4.0*atan(1.0);//pi=PI;

 

    if((n%2)==0)/*如果階數n是偶數*/

    {n2=n/2+1;/**/

     mid=1;//

    }

    else

    {n2=n/2;//n是奇數,則窗口長度為偶數

     mid=0;//

    }

    delay=n/2.0;

    wc1=pi*fln;//

    if(band>=3) wc2=pi*fhn;/*先判斷用戶輸入的數據，如果band參數大於3*/

    switch(band)

    {case 1:

                {for(i=0;i<=n2;i++)

                 {s=i-delay;//

                  h[i]=(sin(wc1*s/fs)/(pi*s))*window(wn,n+1,i,beta);//低通,窗口長度=階數+1，故為n+1

                  h[n-i]=h[i];

                      }

                  if(mid==1) h[n/2]=wc1/pi;//n為偶數時，修正中間值系數

                  break;

                   }

 

    case 2:

                {for(i=0;i<n2;i++)

                {s=i-delay;

                 h[i]=(sin(wc2*s/fs)-sin(wc1*s/fs))/(pi*s);//帶通-//對

                 h[i]=h[i]*window(wn,n+1,i,beta);

                 h[n-i]=h[i];

                   }

                   if(mid==1)h[n/2]=(wc2-wc1)/pi;//對

                  break;

                   }

    case 3:

                   {for(i=0;i<=n2;i++)

                   {s=i-delay;

                   h[i]=(sin(wc1*s/fs)+sin(pi*s)-sin(wc2*s/fs))/(pi*s);//帶阻-//對

                   h[i]=h[i]*window(wn,n+1,i,beta);

                   h[n-i]=h[i];

                    }

                   if(mid==1)h[n/2]=(wc1+pi-wc2)/pi;

                   break;

                   }

    case 4:    

             { for(i=0;i<=n2;i++)

               {s=i-delay;

                h[i]=(sin(pi*s)-sin(wc1*s/fs))/(pi*s);//高通-//對

                h[i]=h[i]*window(wn,n+1,i,beta);

                 h[n-i]=h[i];

                }

               if(mid==1) h[n/2]=1.0-wc1/pi;//對

               break;

              }

                     }

//    for (int _i=0;_i<n+1;_i++)

//    {

//        h[_i]*=(double)(n+1);

//    }

}

double CFIRWIN::window(int type,int n,int i,double beta)

{

    int k;

    double pi,w;

    pi=4.0*atan(1.0);//pi=PI;

    w=1.0;

 

    switch(type)

    {case 1:

    {w=1.0;//矩形窗

    break;

    }

    case 2:

    {k=(n-2)/10;

    if(i<=k)

    w=0.5*(1.0-cos(i*pi/(k+1)));//圖基窗

    break;

    }

    case 3:

    {w=1.0-fabs(1.0-2*i/(n-1.0));//三角窗

    break;

    }

    case 4:

    {w=0.5*(1.0-cos(2*i*pi/(n-1)));//漢寧窗

    break;

    }

    case 5:

    {w=0.54-0.46*cos(2*i*pi/(n-1));//海明窗

    break;

    }

    case 6:

    {w=0.42-0.5*cos(2*i*pi/(n-1))+0.08*cos(4*i*pi/(n-1));//布萊克曼窗

    break;

    }

    case 7:

    {w=kaiser(i,n,beta);//凱塞窗

    break;

    }

    }

    return(w);

}

double CFIRWIN:: kaiser(int i,int n,double beta)

{

    double a,w,a2,b1,b2,beta1;

    

    b1=bessel0(beta);

    a=2.0*i/(double)(n-1)-1.0;

    a2=a*a;

    beta1=beta*sqrt(1.0-a2);

    b2=bessel0(beta1);

    w=b2/b1;

    return(w);

}    

double CFIRWIN::bessel0(double x)

{

    int i;

    double d,y,d2,sum;

    y=x/2.0;

    d=1.0;

    sum=1.0;

    for(i=1;i<=25;i++)

    {d=d*y/i;

    d2=d*d;

    sum=sum+d2;

    if(d2<sum*(1.0e-8)) break;

    }

    return(sum);

}

 

利用firwin這個函數，我們就可以得到hd(n)了。接下來的工作就是對輸入數據序列進行濾波了，由第一部分的公式可以知道，此時有兩種做法。
1.直接按照卷積公式進行計算
2.利用FFT先將x(n)和hd(n)變換到頻域上，得到X(K)和H(k)后相乘得到Y(K)，再進行IFFT即可得到y(n)
下面給出具體代碼：

void CWaveProcess::Filter(float *pfSignal,DWORD dwLenSignal,double *h,int N)

{     

    //法1，直接計算卷積

    double *Input_Buffer;

    double Output_Data = 0;

    Input_Buffer = (double *) malloc(sizeof(double)*N);  

    memset(Input_Buffer,

           0,

           sizeof(double)*N);

    int Count = 0;

    while(1)

    {   

        if(Count==dwLenSignal) break;

         

        Save_Input_Date (pfSignal[Count],

                     N,

                     Input_Buffer);

       

        Output_Data = Real_Time_FIR_Filter(h,

                                           N,

                                           Input_Buffer);

        pfSignal[Count]=Output_Data;

        Count++;

    }

    //法2，傅里葉變換相乘后，做反傅里葉變換

/*    int nPower=(int)(log(N)/log(2))+1;

    int nLen=1<<nPower;

    Complex *A=new Complex[nLen];

    Complex *B=new Complex[nLen];

    Complex *C=new Complex[nLen];

    int nBlock = (dwLenSignal+nLen-1)/nLen;

    CFFT *pA=new CFFT;

    CFFT *pB=new CFFT;

    CFFT *pC=new CFFT;

    for(int i=0; i<nBlock; i++)

    {

        for(int j=0; j<nLen; j++)

        {

            if ((DWORD)(i*nLen+j)<dwLenSignal)

            {

                A[j].real=pfSignal[(DWORD)(i*nLen+j)];

                A[j].imag=0.0;

            }

            else 

            {



                A[j].real=0.0;

                A[j].imag=0.0;

            }

            if (j<N)

            {

                B[j].real=h[j];

                B[j].imag=0.0;

            }

            else 

            {

                B[j].real=0.0;

                B[j].imag=0.0;

            }

        }

        pA->MYFFT(A,nLen);

        pB->MYFFT(B,nLen);

        for(int _i=0;_i<nLen;_i++)

        {

            C[_i]=A[_i]*B[_i];//在頻域進行乘積

        }

        pC->MYFFT(C,nLen,true);//然后再在頻域反變換回時域，就是卷積

        for(j=0;j<nLen;j++)

        {

            if ((DWORD)(i*nLen+j)<dwLenSignal)

                pfSignal[(DWORD)(i*nLen+j)]=C[j].real;

            else break;

        }    

    }

    

    

    delete pA;

    delete pB;

    delete pC;*/

}

double CWaveProcess::Real_Time_FIR_Filter(double *b,

                            int     b_Lenth,

                            double *Input_Data)

{    

    int Count;

    double Output_Data = 0;

    

    Input_Data += b_Lenth - 1;  

    

   for(Count = 0; Count < b_Lenth ;Count++)

    {        

            Output_Data += (*(b + Count)) *

                            (*(Input_Data - Count));                         

    }    

    return (double)Output_Data;

}

void CWaveProcess::Save_Input_Date (double Scand,

                      int    Depth,

                      double *Input_Data)

{

    int Count;

  

    for(Count = 0 ; Count < Depth-1 ; Count++)

    {

        *(Input_Data + Count) = *(Input_Data + Count + 1);

    }

    

    *(Input_Data + Depth-1) = Scand;

}

 

 

實際對比兩種方法，發現通過fft算法來濾波可提高速度。下面貼出fft算法實現過程，基本思路是，逆序，蝶形計算，利用三重循環控制實現。

void CFFT::MYFFT(Complex *A, int N, bool ifft)//當給ifft賦真值的話進行反變換

{

    Complex T;

    int m=(int)(log(N)/log(2)),k,P,B,j=N/2,L,i;//m是級數，為2為底N的對數

    

    for(i=1;i<=N-2;i++)//倒序實現，因為經過m次偶奇抽選之后，先前順序被打亂了，但是打亂后的順序是有規律的

    {    

        if(i<j)

        {

            T=A[i];

            A[i]=A[j];

            A[j]=T;

        }

        k=N/2;

        while(j>=k)

        {

            j-=k;

            k=k/2;

        }

        j+=k;

    }

    

    for(L=1;L<=m;L++)//FFT實現（核心算法時域抽取法）

    {    

        B=1<<(L-1);//這是2的L-1次方

        for(j=0;j<=B-1;j++)

        {

            P=(1<<(m-L))*j;

            if(ifft==false)//默認算法為傅里葉正變換

            {

                for(k=j;k<=N-1;k+=1<<L)

                {

                    T=A[k]+A[k+B]*Complex(cos(-2*PI*P/N),sin(-2*PI*P/N));

                    A[k+B]=A[k]-A[k+B]*Complex(cos(-2*PI*P/N),sin(-2*PI*P/N));

                    A[k]=T;

                }

            }

            else//ifft為真的話進行傅里葉反變換

            {

                for(k=j;k<=N-1;k+=1<<L)

                {

                    T=A[k]+A[k+B]*Complex(cos(-2*PI*P/N),sin(2*PI*P/N));//反變換得取共軛

                    A[k+B]=A[k]-A[k+B]*Complex(cos(-2*PI*P/N),sin(2*PI*P/N));

                    A[k]=T;

                }

            }

        }

    }

    if(ifft==true)//反變換還得除以N

    {

        for(k=0;k<N;k++)

            A[k]=A[k]/N;

    }

}

 

3.結束語

進行數字信號處理可利用的工具有很多，比如matlab，LabVIEW等，這些工具都很強大，使用起來也特別方便。通常C語言要用大量代碼實現的過程，matlab一句代碼，LabVIEW一個圖形就可以代替，因為已經做好了封裝，方便使用。但是用C語言的好處就是，能對底層進行修改，使程序設計更加靈活。同時，進行底層語言的編寫，可以深入理解原理，加深對數字信號處理這門課程基礎知識的掌握。
————————————————
版權聲明：本文為CSDN博主「hitwhzhongqiu」的原創文章，遵循 CC 4.0 BY-SA 版權協議，轉載請附上原文出處鏈接及本聲明。
原文鏈接：https://blog.csdn.net/hitwhzhongqiu/article/details/41948623