Implementación del filtro paso banda Butterworth en C++

Bienvenido a proyecto on line, ahora vas a hallar la resolución a lo que necesitas.

Solución:

Sé que esta es una publicación en un hilo antiguo, y normalmente dejaría esto como un comentario, pero aparentemente no puedo hacerlo.

En cualquier caso, para las personas que buscan un 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 algún 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 hay un fragmento de código que ya calcula dicho factor de escala para usted.

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

Finalmente lo encontré. Solo necesito implementar el siguiente código desde el código fuente de matlab a c++. "the_mandrill" tenía razón, necesito agregar la constante de normalización en el coeficiente:

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

EDITAR:
y aquí está la edición final, en la que todo el código devolverá números exactamente iguales a MATLAB:

double *ComputeNumCoeffs(int FilterOrder,double Lcutoff, double Ucutoff, double *DenC)

    double *TCoeffs;
    double *NumCoeffs;
    std::complex *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 *)calloc( 2*FilterOrder+1, sizeof(std::complex) );
    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 result = std::complex(-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;

Comentarios y valoraciones

Al final de la página puedes encontrar las notas de otros sys admins, tú todavía tienes el poder dejar el tuyo si te apetece.

---

---

Utiliza Nuestro Buscador