pasmowego realizacja Filtr Butterwortha w C++

Question

ja wdrożenie algorytm analizy obrazu za pomocą OpenCV i C++, ale okazało się, OpenCV robi się żadnej funkcji dla Butterworth pasmowego filtrować oficjalnie. w moim projekcie muszę przejść szereg czasowy pikseli do filtra Butterwortha rzędu 5 i funkcja zwróci przefiltrowane serii pikseli czasowych. Butterworth (pixelseries, kolejność, częstotliwość), jeśli masz pomysł, aby mi pomóc, jak zacząć, daj mi znać. Dziękuję

EDIT: po uzyskaniu pomocy, w końcu wymyślić poniższego kodu. które mogą obliczyć współczynniki licznik i mianownik współczynników ale problemem jest to, że niektóre z tych liczb nie jest tak samo jak wyniki Matlab. tutaj jest mój kod:

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

uzyskać to resutls dla mojego kodu:

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

jednak jeśli chce się testować CO efficinets w tym samym paśmie częstotliwości w programie MATLAB, mam następujące wyniki:

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

mam testu tej stronie: http: //www.exstrom .com/journal/sigproc/code, ale wynik jest taki sam jak mój, nie matlab. ktokolwiek wie dlaczego? lub w jaki sposób mogę uzyskać ten sam wynik, co przybornik matlab?

Answer 1

końcu znalazłem go. po prostu trzeba zaimplementować następujący kod z kodu źródłowego Matlab do C++. „The_mandrill” mieli rację, muszę dodać normalizujące stałej do współczynnika:

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

EDIT: i tu jest ostateczna wersja, która cały kod powróci numery dokładnie równa 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

Wiem, że jest to post na starym wątku, a ja zwykle zostawiam to jako komentarz, ale najwyraźniej nie jestem w stanie tego zrobić.

W każdym przypadku, dla osób szukających podobnego kodu, pomyślałem, że zamieściłbym łącze z tego, skąd pochodzi ten kod (ma również kod C dla innych typów współczynników filtru Butterwortha i jakiś inny kod przetwarzania chłodnego sygnału).

kod znajduje się tutaj: http://www.exstrom.com/journal/sigproc/

Dodatkowo, myślę, że jest to fragment kodu, który oblicza powiedział współczynnik skalowania dla Ciebie już.

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