paso de banda implementación del filtro butterworth en C++

Question

Estoy implementando un algoritmo de análisis de imágenes usando openCV y C++, pero descubrí que openCV no tiene ninguna función para el filtro Butterworth Bandpass oficialmente. en mi proyecto Tengo que pasar una serie temporal de píxeles en el filtro de orden Butterworth 5 y la función devolverá los píxeles de la serie temporal filtrada. Butterworth (pixelseries, orden, frecuencia), si tiene alguna idea para ayudarme a comenzar, háganmelo saber. Gracias

EDIT: Después de recibir ayuda, finalmente se me ocurre el siguiente código. que puede calcular los Coeficientes del Numerador y los Coeficientes del Denominador, pero el problema es que algunos de los números no son tan iguales como los resultados del Matlab. aquí está mi código:

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

consigo esto los resultados de mi código:

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

pero si quiero para probar la co efficinets en misma banda de frecuencia en MATLAB, que obtienen los resultados siguientes:

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

que tienen prueba de este sitio web: http: //www.exstrom .com/journal/sigproc/code, pero el resultado es igual que el mío, no matlab. alguien sabe por qué? o ¿cómo puedo obtener el mismo resultado que la caja de herramientas de matlab?

Answer 1

Finalmente lo encontré. Solo necesito implementar el siguiente código desde el código fuente de matlab a C++. "The_mandrill" tenían razón, tengo que añadir la constante de normalización en el coeficiente:

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

EDIT: y aquí está la edición final, que todo el código volverá números exactamente igual a 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

Sé que esta es una publicación en un hilo viejo, y generalmente dejaría esto como un comentario, pero aparentemente no puedo hacer eso. En cualquier caso, para las personas que buscan código similar, pensé en publicar el enlace desde donde se origina este código (también tiene código C para otros tipos de coeficientes de filtro Butterworth y otro código de procesamiento de señal interesante).

El código se encuentra aquí: http://www.exstrom.com/journal/sigproc/

Además, creo que es una pieza de código que calcula dicho factor de escala para los que ya.

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