在C++中實現帶通butterworth過濾器

Question

我正在實現一個使用openCV和C++的圖像分析算法，但是我發現openCV並沒有正式的Butterworth Bandpass過濾器的任何功能。 在我的項目中，我必須將時間序列的像素傳遞到巴特沃斯5階濾波器，並且該函數將返回濾波後的時間序列像素。巴特沃斯（像素系列，訂單，頻率），如果您有任何想法來幫助我如何開始請讓我知道。謝謝

編輯： 得到幫助後，最後我拿出下面的代碼。它可以計算分子係數和分母系數，但問題是有些數字與matlab結果不一樣。這裏是我的代碼：

#include <iostream> 
#include <stdio.h> 
#include <vector> 
#include <math.h> 

using namespace std; 

#define N 10 //The number of images which construct a time series for each pixel 
#define PI 3.14159 

double *ComputeLP(int FilterOrder) 
{ 
    double *NumCoeffs; 
    int m; 
    int i; 

    NumCoeffs = (double *)calloc(FilterOrder+1, sizeof(double)); 
    if(NumCoeffs == NULL) return(NULL); 

    NumCoeffs[0] = 1; 
    NumCoeffs[1] = FilterOrder; 
    m = FilterOrder/2; 
    for(i=2; i <= m; ++i) 
    { 
     NumCoeffs[i] =(double) (FilterOrder-i+1)*NumCoeffs[i-1]/i; 
     NumCoeffs[FilterOrder-i]= NumCoeffs[i]; 
    } 
    NumCoeffs[FilterOrder-1] = FilterOrder; 
    NumCoeffs[FilterOrder] = 1; 

    return NumCoeffs; 
} 

double *ComputeHP(int FilterOrder) 
{ 
    double *NumCoeffs; 
    int i; 

    NumCoeffs = ComputeLP(FilterOrder); 
    if(NumCoeffs == NULL) return(NULL); 

    for(i = 0; i <= FilterOrder; ++i) 
     if(i % 2) NumCoeffs[i] = -NumCoeffs[i]; 

    return NumCoeffs; 
} 

double *TrinomialMultiply(int FilterOrder, double *b, double *c) 
{ 
    int i, j; 
    double *RetVal; 

    RetVal = (double *)calloc(4 * FilterOrder, sizeof(double)); 
    if(RetVal == NULL) return(NULL); 

    RetVal[2] = c[0]; 
    RetVal[3] = c[1]; 
    RetVal[0] = b[0]; 
    RetVal[1] = b[1]; 

    for(i = 1; i < FilterOrder; ++i) 
    { 
     RetVal[2*(2*i+1)] += c[2*i] * RetVal[2*(2*i-1)] - c[2*i+1] * RetVal[2*(2*i-1)+1]; 
     RetVal[2*(2*i+1)+1] += c[2*i] * RetVal[2*(2*i-1)+1] + c[2*i+1] * RetVal[2*(2*i-1)]; 

     for(j = 2*i; j > 1; --j) 
     { 
      RetVal[2*j] += b[2*i] * RetVal[2*(j-1)] - b[2*i+1] * RetVal[2*(j-1)+1] + 
       c[2*i] * RetVal[2*(j-2)] - c[2*i+1] * RetVal[2*(j-2)+1]; 
      RetVal[2*j+1] += b[2*i] * RetVal[2*(j-1)+1] + b[2*i+1] * RetVal[2*(j-1)] + 
       c[2*i] * RetVal[2*(j-2)+1] + c[2*i+1] * RetVal[2*(j-2)]; 
     } 

     RetVal[2] += b[2*i] * RetVal[0] - b[2*i+1] * RetVal[1] + c[2*i]; 
     RetVal[3] += b[2*i] * RetVal[1] + b[2*i+1] * RetVal[0] + c[2*i+1]; 
     RetVal[0] += b[2*i]; 
     RetVal[1] += b[2*i+1]; 
    } 

    return RetVal; 
} 

double *ComputeNumCoeffs(int FilterOrder) 
{ 
    double *TCoeffs; 
    double *NumCoeffs; 
    int i; 

    NumCoeffs = (double *)calloc(2*FilterOrder+1, sizeof(double)); 
    if(NumCoeffs == NULL) return(NULL); 

    TCoeffs = ComputeHP(FilterOrder); 
    if(TCoeffs == NULL) return(NULL); 

    for(i = 0; i < FilterOrder; ++i) 
    { 
     NumCoeffs[2*i] = TCoeffs[i]; 
     NumCoeffs[2*i+1] = 0.0; 
    } 
    NumCoeffs[2*FilterOrder] = TCoeffs[FilterOrder]; 

    free(TCoeffs); 

    return NumCoeffs; 
} 

double *ComputeDenCoeffs(int FilterOrder, double Lcutoff, double Ucutoff) 
{ 
    int k;   // loop variables 
    double theta;  // PI * (Ucutoff - Lcutoff)/2.0 
    double cp;  // cosine of phi 
    double st;  // sine of theta 
    double ct;  // cosine of theta 
    double s2t;  // sine of 2*theta 
    double c2t;  // cosine 0f 2*theta 
    double *RCoeffs;  // z^-2 coefficients 
    double *TCoeffs;  // z^-1 coefficients 
    double *DenomCoeffs;  // dk coefficients 
    double PoleAngle;  // pole angle 
    double SinPoleAngle;  // sine of pole angle 
    double CosPoleAngle;  // cosine of pole angle 
    double a;   // workspace variables 

    cp = cos(PI * (Ucutoff + Lcutoff)/2.0); 
    theta = PI * (Ucutoff - Lcutoff)/2.0; 
    st = sin(theta); 
    ct = cos(theta); 
    s2t = 2.0*st*ct;  // sine of 2*theta 
    c2t = 2.0*ct*ct - 1.0; // cosine of 2*theta 

    RCoeffs = (double *)calloc(2 * FilterOrder, sizeof(double)); 
    TCoeffs = (double *)calloc(2 * FilterOrder, sizeof(double)); 

    for(k = 0; k < FilterOrder; ++k) 
    { 
     PoleAngle = PI * (double)(2*k+1)/(double)(2*FilterOrder); 
     SinPoleAngle = sin(PoleAngle); 
     CosPoleAngle = cos(PoleAngle); 
     a = 1.0 + s2t*SinPoleAngle; 
     RCoeffs[2*k] = c2t/a; 
     RCoeffs[2*k+1] = s2t*CosPoleAngle/a; 
     TCoeffs[2*k] = -2.0*cp*(ct+st*SinPoleAngle)/a; 
     TCoeffs[2*k+1] = -2.0*cp*st*CosPoleAngle/a; 
    } 

    DenomCoeffs = TrinomialMultiply(FilterOrder, TCoeffs, RCoeffs); 
    free(TCoeffs); 
    free(RCoeffs); 

    DenomCoeffs[1] = DenomCoeffs[0]; 
    DenomCoeffs[0] = 1.0; 
    for(k = 3; k <= 2*FilterOrder; ++k) 
     DenomCoeffs[k] = DenomCoeffs[2*k-2]; 


    return DenomCoeffs; 
} 

void filter(int ord, double *a, double *b, int np, double *x, double *y) 
{ 
    int i,j; 
    y[0]=b[0] * x[0]; 
    for (i=1;i<ord+1;i++) 
    { 
     y[i]=0.0; 
     for (j=0;j<i+1;j++) 
      y[i]=y[i]+b[j]*x[i-j]; 
     for (j=0;j<i;j++) 
      y[i]=y[i]-a[j+1]*y[i-j-1]; 
    } 
    for (i=ord+1;i<np+1;i++) 
    { 
     y[i]=0.0; 
     for (j=0;j<ord+1;j++) 
      y[i]=y[i]+b[j]*x[i-j]; 
     for (j=0;j<ord;j++) 
      y[i]=y[i]-a[j+1]*y[i-j-1]; 
    } 
} 




int main(int argc, char *argv[]) 
{ 
    //Frequency bands is a vector of values - Lower Frequency Band and Higher Frequency Band 

    //First value is lower cutoff and second value is higher cutoff 
    double FrequencyBands[2] = {0.25,0.375};//these values are as a ratio of f/fs, where fs is sampling rate, and f is cutoff frequency 
    //and therefore should lie in the range [0 1] 
    //Filter Order 

    int FiltOrd = 5; 

    //Pixel Time Series 
    /*int PixelTimeSeries[N]; 
    int outputSeries[N]; 
    */ 
    //Create the variables for the numerator and denominator coefficients 
    double *DenC = 0; 
    double *NumC = 0; 
    //Pass Numerator Coefficients and Denominator Coefficients arrays into function, will return the same 

    NumC = ComputeNumCoeffs(FiltOrd); 
    for(int k = 0; k<11; k++) 
    { 
     printf("NumC is: %lf\n", NumC[k]); 
    } 
    //is A in matlab function and the numbers are correct 
    DenC = ComputeDenCoeffs(FiltOrd, FrequencyBands[0], FrequencyBands[1]); 
    for(int k = 0; k<11; k++) 
    { 
     printf("DenC is: %lf\n", DenC[k]); 
    } 
    double y[5]; 
    double x[5]={1,2,3,4,5}; 
    filter(5, DenC, NumC, 5, x, y);  
    return 1; 
} 

我得到這個resutls我的代碼：

B = 1,0，-5,0,10,0，-10,0,5,0， -1 A = 1.000000000000000，-4.945988709743181，13.556489496973796，-24.700711850327743， 32.994881546824828，-33.180726698160655，25.546126213403539，-14.802008410165968， 6.285430089797051，-1.772929809750849，0.277753012228403

，但如果我想測試合作在MATLAB相同頻帶efficinets時，得到以下結果：

>> [B, A]=butter(5, [0.25,0.375]) 

B = 0.0002，0，-0.0008，0，0.0016，0，-0.0016，0，0.0008，0，-0.0002

A = 1.0000，-4.9460，13.5565，-24.7007，32.9948，-33.1806，25.5461，-14.8020，6.2854，-1.7729，0.2778

我有測試此網址：http：//www.exstrom .com/journal/sigproc/code，但結果與我的相同，不是matlab。有人知道爲什麼嗎？或者我如何能得到與matlab工具箱相同的結果？

Answer 1

我終於找到它了。 我只需要將matlab源代碼中的以下代碼實現爲C++。 「the_mandrill」是對的，我需要的標準化常數添加到係數：

kern = exp(-j*w*(0:length(b)-1)); 
b = real(b*(kern*den(:))/(kern*b(:))); 

編輯： 這裏是最終版本，這整個代碼將完全等於回到號碼MATLAB：

double *ComputeNumCoeffs(int FilterOrder,double Lcutoff, double Ucutoff, double *DenC) 
{ 
    double *TCoeffs; 
    double *NumCoeffs; 
    std::complex<double> *NormalizedKernel; 
    double Numbers[11]={0,1,2,3,4,5,6,7,8,9,10}; 
    int i; 

    NumCoeffs = (double *)calloc(2*FilterOrder+1, sizeof(double)); 
    if(NumCoeffs == NULL) return(NULL); 

    NormalizedKernel = (std::complex<double> *)calloc(2*FilterOrder+1, sizeof(std::complex<double>)); 
    if(NormalizedKernel == NULL) return(NULL); 

    TCoeffs = ComputeHP(FilterOrder); 
    if(TCoeffs == NULL) return(NULL); 

    for(i = 0; i < FilterOrder; ++i) 
    { 
     NumCoeffs[2*i] = TCoeffs[i]; 
     NumCoeffs[2*i+1] = 0.0; 
    } 
    NumCoeffs[2*FilterOrder] = TCoeffs[FilterOrder]; 
    double cp[2]; 
    double Bw, Wn; 
    cp[0] = 2*2.0*tan(PI * Lcutoff/ 2.0); 
    cp[1] = 2*2.0*tan(PI * Ucutoff/2.0); 

    Bw = cp[1] - cp[0]; 
    //center frequency 
    Wn = sqrt(cp[0]*cp[1]); 
    Wn = 2*atan2(Wn,4); 
    double kern; 
    const std::complex<double> result = std::complex<double>(-1,0); 

    for(int k = 0; k<11; k++) 
    { 
     NormalizedKernel[k] = std::exp(-sqrt(result)*Wn*Numbers[k]); 
    } 
    double b=0; 
    double den=0; 
    for(int d = 0; d<11; d++) 
    { 
     b+=real(NormalizedKernel[d]*NumCoeffs[d]); 
     den+=real(NormalizedKernel[d]*DenC[d]); 
    } 
    for(int c = 0; c<11; c++) 
    { 
     NumCoeffs[c]=(NumCoeffs[c]*den)/b; 
    } 

    free(TCoeffs); 
    return NumCoeffs; 
} 

Answer 2

我知道這是一箇舊線程的帖子，我通常會將此留作評論，但我顯然無法做到這一點。在任何情況下，對於尋找類似代碼的人來說，我認爲我會發布此代碼來自何處的鏈接（它還具有用於其他類型Butterworth濾波器係數和一些其他酷信號處理代碼的C代碼）。

的代碼位於： http://www.exstrom.com/journal/sigproc/

此外，我認爲這是一段代碼，它計算說已經爲您縮放因子。

/********************************************************************** 
sf_bwbp - calculates the scaling factor for a butterworth bandpass filter. 
The scaling factor is what the c coefficients must be multiplied by so 
that the filter response has a maximum value of 1. 
*/ 

double sf_bwbp(int n, double f1f, double f2f) 
{ 
    int k;   // loop variables 
    double ctt;  // cotangent of theta 
    double sfr, sfi; // real and imaginary parts of the scaling factor 
    double parg;  // pole angle 
    double sparg;  // sine of pole angle 
    double cparg;  // cosine of pole angle 
    double a, b, c; // workspace variables 

    ctt = 1.0/tan(M_PI * (f2f - f1f)/2.0); 
    sfr = 1.0; 
    sfi = 0.0; 

    for(k = 0; k < n; ++k) 
    { 
     parg = M_PI * (double)(2*k+1)/(double)(2*n); 
     sparg = ctt + sin(parg); 
     cparg = cos(parg); 
     a = (sfr + sfi)*(sparg - cparg); 
     b = sfr * sparg; 
     c = -sfi * cparg; 
     sfr = b - c; 
     sfi = a - b - c; 
    } 

    return(1.0/sfr); 
}