¿Cómo funciona un filtro de paso bajo mediante programación?

He estado trabajando en un filtro de paso bajo simple para la medición de <100 Hz en mi aplicación. Pero hasta ahora, estoy luchando con la teoría detrás de todo. Es genial que lo haya hecho funcionar, pero realmente lo disfrutaría si supiera cómo / por qué está funcionando.

Encontré el siguiente código:

void getLPCoefficientsButterworth2Pole(const int samplerate, const double cutoff, double* const ax, double* const by)
{
    double PI = M_PI;
    double sqrt2 = sqrt(2);

    double QcRaw  = (2 * PI * cutoff) / samplerate; // Find cutoff frequency in [0..PI]
    double QcWarp = tan(QcRaw); // Warp cutoff frequency

    double gain = 1 / ( 1 + sqrt2 / QcWarp + 2 / ( QcWarp * QcWarp ) );

    by[2] = ( 1 - sqrt2 / QcWarp + 2 / ( QcWarp * QcWarp ) ) * gain;
    by[1] = ( 2 - 2 * 2 / ( QcWarp * QcWarp ) ) * gain;
    by[0] = 1;

    ax[0] = 1 * gain;
    ax[1] = 2 * gain;
    ax[2] = 1 * gain;
}

Para calcular los coeficientes. Luego, en las muestras de audio, las 'paso bajo' de esta manera:

        xv[2] = xv[1];
        xv[1] = xv[0];

        xv[0] = pData[j];
        yv[2] = yv[1];
        yv[1] = yv[0];

        yv[0] = (ax[0] * xv[0] + ax[1] * xv[1] + ax[2] * xv[2]
                   - by[1] * yv[0]
                   - by[2] * yv[1]);

        pData[j] = yv[0];

Para obtener un diseño de paso bajo.

Me pregunto algunas cosas:

1. Recibo las muestras de audio en una simple matriz flotante *. ¿Cuál es ese número flotante? Lo único que veo es un número, ¿cómo es ese sonido?
2. El código está usando cálculos pasados ​​(tres de ellos) en el nuevo cálculo por muestra. ¿Eso significa que las primeras 2 muestras de datos no se filtran correctamente? (No es que importe porque son solo 2 muestras, pero me pregunto)
3. Intentando aprender todo, encontré un par de fórmulas para el filtro Butterworth (2º Polo). ¿Cómo se reflejan esas fórmulas en este código? Ninguna de las fórmulas que encontré tiene estos cálculos que puede ver en la función 'getLPCoefficientsButterworth2Pole ()'.

1. El número flotante * matriz es un puntero a la matriz. Es un número único que contiene la dirección del primer elemento de la matriz de valores flotantes.

2. Por lo general, la condición inicial (es decir, los elementos 'pasados' iniciales de x e y) son 0, pero si sus valores no son iguales a 0 tampoco es un gran problema, porque después de un tiempo las condiciones iniciales no tienen efecto en la salida señal para cualquier filtro estable. Y su filtro es obviamente estable.

3. Una función de transferencia de paso bajo de segundo orden con la característica de Butterworth y la frecuencia de corte (con en hercios) viene dada por$\omega_a=2*\pi*f_a$$f_a$

$$H(s) = \frac{\omega_a^2}{s^2+\sqrt{2}\omega_a s + \omega_a^2}\tag{1}$$ Este es un resultado estándar que puede encontrar fácilmente en la web . Para obtener un filtro de tiempo discreto, puede aplicar una llamada transformación bilineal:

$$s = 2f_s\frac{z-1}{z+1}\tag{2}$$

donde es la frecuencia de muestreo. Esto es necesario porque el eje de frecuencia de las señales analógicas ( ) debe asignarse al rango de frecuencia permitido para señales de tiempo discreto ( ). Dado que esta transformación deforma las frecuencias, necesitamos calcular la frecuencia de corte analógica deseada a partir de la frecuencia de corte dada en el dominio de tiempo discreto:$f_s$$0\le f\le\infty$$0\le f\le f_s/2$$f_d$

$$\omega_a = 2f_s \tan(\frac{\omega_d}{2})\text{ with }\omega_d = 2\pi f_d/f_s$$

Si inserta (2) en (1) obtendrá

$$H(z) = k\frac{z^2+b_1z+b_2}{z^2+a_1z+a_2}\tag{3}$$

con

$$k=\frac{\alpha^2}{1+\sqrt{2}\alpha+\alpha^2}\text{ (gain)}$$ $$a_1 = \frac{2(\alpha^2-1)}{1+\sqrt{2}\alpha+\alpha^2},\quad a_2=\frac{1-\sqrt{2}\alpha+\alpha^2}{1+\sqrt{2}\alpha+\alpha^2},\quad b_1=2,\quad b_2=1$$

donde he usado . La función de transferencia (3) corresponde a la ecuación del dominio del tiempo$\alpha=\tan(\frac{\omega_d}{2})$

$$y(n) = kx(n) + kb_1x(n-1) + kb_2x(n-2) - a_1y(n-1) - a_2y(n-2)$$

Como puede ver, hay bastante parecido con su código. Sin embargo, también hay algunas diferencias y debe verificar la respuesta de frecuencia de su filtro. Creo que las ecuaciones anteriores son correctas, pero depende de usted verificarlas y decidir qué versión es la correcta.

Preguntó cómo funciona un filtro de paso bajo y mencionó que el filtro utiliza valores pasados ​​de sus datos. Esta es una discusión no técnica de lo que sucede en un filtro de paso bajo.

El filtro de paso bajo toma diferentes vistas (desplazadas en el tiempo) de su señal, las escala y las agrega. Puede imaginarse dibujando su señal 3 veces, una siendo actual, la segunda desplazada por un tiempo de muestra, la tercera desplazada por 2 veces de muestra.

A bajas frecuencias, todas las vistas se ven muy similares (el cambio por una sola muestra apenas cambia la ubicación de la señal en cualquier momento). En este caso, las tres versiones se sumarán de manera constructiva (o al menos no destructiva), de modo que la señal pase a través del filtro.

Ahora moviéndose a frecuencias más altas, cada versión desplazada de la señal se vuelve más distinta en cualquier instante dado (punto de muestra) y, de hecho, incluso puede invertir el signo. A estas frecuencias más altas, las tres versiones de su señal tienden a cancelarse (agregadas destructivamente) por lo que la señal se atenúa.

Los diferentes tipos de filtros organizan este tipo de interferencia constructiva / destructiva en bandas de frecuencia apropiadas para crear filtros de paso bajo, paso de banda o paso alto.

Depende de cómo lo quieras. Para mí lo he implementado en C. Entonces, lo que hice fue que generé coeficientes de filtro en el matlabuso firls. Luego almacené esos coeficientes de filtro en una matriz y luego pasé estos coeficientes de filtro a la función de filtro en C. En matlabél es bastante fácil generar los coeficientes del filtro, y luego usar la operación de convolución o usar un filtro específico. Pero también Cdebes implementar el filtro. Por mi parte, tengo que hacerlo para DSP, así que lo implementé Cy obtuve los resultados al matlabusar la función mex. Para calcular los coeficientes, use este comando en matlab:

n=17                  %filter order
f=[0 0.167 0.333 1]   %Frequency band edges
a=[1 1 0 0]           %Desired amplitudes
fir= firls (n, f,a )