thực thi bộ lọc dải băng trong C++

Question

Tôi đang triển khai thuật toán phân tích hình ảnh bằng cách sử dụng openCV và C++, nhưng tôi phát hiện ra openCV không có bất kỳ chức năng nào cho bộ lọc Butterworth Bandpass chính thức. trong dự án của tôi, tôi phải chuyển một chuỗi thời gian pixel vào bộ lọc thứ tự Butterworth 5 và hàm sẽ trả về các pixel chuỗi thời gian đã lọc. Butterworth (pixeleries, thứ tự, tần số), nếu bạn có bất kỳ ý tưởng nào để giúp tôi biết cách bắt đầu, hãy cho tôi biết. Cảm ơn bạn

EDIT: sau khi nhận trợ giúp, cuối cùng tôi tìm ra mã sau đây. có thể tính toán hệ số Numerator và hệ số mẫu số, nhưng vấn đề là một số các số không giống như kết quả MATLAB. đây là mã của tôi:

#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; 
} 

tôi có được điều này resutls cho mã của tôi:

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

nhưng nếu tôi muốn thử nghiệm đồng efficinets trong ban nhạc cùng một tần số trong MATLAB, tôi nhận được kết quả sau:

>> [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

tôi có kiểm tra trang web này: http: //www.exstrom .com/journal/sigproc/code, nhưng kết quả bằng như của tôi, không phải MATLAB. ai biết tại sao? hoặc làm thế nào tôi có thể nhận được kết quả tương tự như hộp công cụ MATLAB?

Answer 1

Cuối cùng tôi đã tìm thấy nó. Tôi chỉ cần triển khai đoạn mã sau từ mã nguồn MATLAB tới C++. "The_mandrill" đã đúng, tôi cần thêm hằng số bình thường hóa vào hệ số:

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

EDIT: và đây là phiên bản cuối cùng, toàn bộ mã sẽ trở lại con số chính xác bằng 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

Tôi biết đây là một bài đăng trên một chủ đề cũ và tôi thường để lại nhận xét này, nhưng tôi dường như không thể làm điều đó.

Trong mọi trường hợp, đối với những người tìm kiếm mã tương tự, tôi nghĩ tôi sẽ đăng liên kết từ nơi mã này bắt nguồn (mã này cũng có mã C cho các loại hệ số lọc Butterworth khác và một số mã xử lý tín hiệu khác).

Mã này nằm ở đây: http://www.exstrom.com/journal/sigproc/

Ngoài ra, tôi nghĩ rằng đó là một đoạn mã này sẽ tính nói rộng yếu tố cho bạn rồi.

/********************************************************************** 
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); 
}