Como um filtro passa-baixo funciona programaticamente?

Eu tenho trabalhado em um filtro passa-baixo simples para medição <100 Hz em meu aplicativo. Mas até agora, estou lutando com a teoria por trás de tudo. É legal que eu tenha funcionado, mas eu realmente apreciaria se soubesse como / por que está funcionando.

Encontrei o seguinte 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 os coeficientes. Então, nas amostras de áudio, eu as transmito da seguinte maneira:

        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 obter um design lowpass.

Estou pensando em algumas coisas:

Eu recebo as amostras de áudio em uma matriz float * simples. Qual é esse número flutuante? A única coisa que vejo é um número, como isso é uma peça de som?
O código está usando cálculos anteriores (três deles) no novo cálculo por amostra. Isso significa que as 2 primeiras amostras de dados não foram filtradas corretamente? (não que isso importe porque são apenas 2 amostras, mas apenas imaginando)
Tentando aprender tudo, encontrei algumas fórmulas para o filtro Butterworth (2º pólo). Como essas fórmulas são refletidas neste código? Nenhuma das fórmulas que encontrei possui esses cálculos que você pode ver na função 'getLPCoefficientsButterworth2Pole ()'.

lowpass-filter c Niek van der Steen
fonte

Não estou tentando ser desrespeitoso aqui, mas sua afirmação "como isso é uma peça de som?" parece indicar que você não entende os princípios básicos do processamento discreto no tempo. Será muito difícil resolver qualquer problema de DSP sem entender o básico, como o teorema da amostragem, quantização, sistema LTI etc. Eu recomendaria algum tempo para um bom livro. Este é gratuito dspguide.com/pdfbook.htm

Hilmar

Respostas:

O número float * array é um ponteiro para o array. É um número único que contém o endereço do primeiro elemento da matriz de valores flutuantes.
Normalmente, a condição inicial (isto é, os elementos iniciais do passado de xey) são 0, mas se seus valores não forem iguais a 0, também não é um grande problema, porque depois de um tempo as condições iniciais não terão efeito na saída. sinal para qualquer filtro estável. E seu filtro é obviamente estável.
Uma função de transferência lowpass de segunda ordem com característica Butterworth e frequência de corte (com em Hertz) é dada por $\omega_a=2*\pi*f_a$ $f_a$

\begin{matrix} (1) & H (s) = \frac{ω_{a}^{2}}{s^{2} + \sqrt{2} ω_{a} s + ω_{a}^{2}} \end{matrix}

$H(s) = \frac{\omega_a^2}{s^2+\sqrt{2}\omega_a s + \omega_a^2}\tag{1}$ Esse é um resultado padrão que você pode encontrar facilmente na web . Para obter um filtro de tempo discreto, você pode aplicar a chamada transformação bilinear:

\begin{matrix} (2) & s = 2 f_{s} \frac{z - 1}{z + 1} \end{matrix}

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

onde é a frequência de amostragem. Isso é necessário porque o eixo de frequência dos sinais analógicos ( ) precisa ser mapeado para a faixa de frequência permitida para sinais de tempo discreto ( ). Como essa transformação distorce as frequências, precisamos calcular a frequência de corte analógica desejada a partir da frequência de corte dada no domínio de tempo discreto: $f_s$ $0\le f\le\infty$ $0\le f\le f_s/2$ $f_d$

ω_{a} = 2 f_{s} \tan (\frac{ω_{d}}{2}) with ω_{d} = 2 π f_{d} / f_{s}

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

Se você inserir (2) em (1), obtém

\begin{matrix} (3) & H (z) = k \frac{z^{2} + b_{1} z + b_{2}}{z^{2} + a_{1} z + a_{2}} \end{matrix}

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

com

k = \frac{α^{2}}{1 + \sqrt{2} α + α^{2}} (gain)

$k=\frac{\alpha^2}{1+\sqrt{2}\alpha+\alpha^2}\text{ (gain)}$

a_{1} = \frac{2 (α^{2} - 1)}{1 + \sqrt{2} α + α^{2}}, a_{2} = \frac{1 - \sqrt{2} α + α^{2}}{1 + \sqrt{2} α + α^{2}}, b_{1} = 2, b_{2} = 1

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

onde eu usei . A função de transferência (3) corresponde à equação do domínio do tempo $\alpha=\tan(\frac{\omega_d}{2})$

y (n) = k x (n) + k b_{1} x (n - 1) + k b_{2} x (n - 2) - a_{1} y (n - 1) - a_{2} y (n - 2)

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

Como você pode ver, há bastante semelhança com o seu código. No entanto, também existem algumas diferenças e você deve verificar a resposta de frequência do seu filtro. Penso que as equações acima estão corretas, mas cabe a você verificá-las e decidir qual versão é a correta.

Matt L.
fonte

Obrigado. O primeiro, eu entendo. Mas ainda é uma variedade de carros alegóricos. Quero dizer: o que um único valor flutuante representa?

Niek van der Steen

Um único valor flutuante é uma das suas amostras de entrada.

Matt L.

Em resumo: o som é uma variação da pressão do ar, convertida por um transdutor (microfone) em um sinal elétrico cujas variações seguem o mesmo padrão da pressão do ar. A amplitude do sinal elétrico é medida N vezes (= amostrador) a cada segundo para produzir a sequência de números que você vê em sua matriz float *. Existem dezenas de apresentações on-line sobre processamento digital (de áudio), explicando esse processo de maneira mais detalhada.

Pichenettes

Editei minha resposta para abordar a última parte da sua pergunta.

Matt L.

Você está usando x (n - algo), o que é x? Presumo que 'n' é a amostra de entrada? Boa resposta!

Niek van der Steen

Você perguntou como um filtro passa-baixo funciona e mencionou que o filtro usa valores passados dos seus dados. Esta é uma discussão não técnica do que acontece em um filtro passa-baixo.

O filtro passa-baixo obtém visões diferentes (deslocadas no tempo) do seu sinal, as escala e as adiciona. Você pode imaginar desenhar seu sinal três vezes, sendo um atual, o segundo sendo alterado por um tempo de amostra, e o terceiro sendo alterado por 2 vezes.

Em frequências baixas, todas as visualizações parecem muito semelhantes (a troca por uma única amostra quase não muda onde você está no sinal a qualquer instante). Nesse caso, as três versões serão adicionadas de maneira construtiva (ou pelo menos não destrutiva), de modo que o sinal passe pelo filtro.

Agora, passando para frequências mais altas, cada versão alterada do sinal se torna mais distinta a qualquer instante (ponto de amostragem) e, de fato, pode até reverter o sinal. Nessas frequências mais altas, as três versões do seu sinal tendem a ser canceladas (adicionadas destrutivamente), para que o sinal seja atenuado.

Diferentes tipos de filtros organizam para que esse tipo de interferência construtiva / destrutiva ocorra nas bandas de frequência apropriadas para criar filtros passa-baixo, passa-banda ou passa-alto.

user2718
fonte

Depende de como você deseja. Para mim eu implementei no C. Então, o que eu fiz foi que eu gerava coeficientes de filtro matlabusando firls. Em seguida, armazenado que os coeficientes de filtro de uma matriz e, em seguida, passar esses coeficientes de filtro para o filtro em função C. Em matlabque é bastante fácil gerado os coeficientes do filtro, e depois usar a operação de convolução ou usar um filtro específico. Mas Cvocê também precisa implementar o filtro. Da minha parte, eu tenho que fazer isso para DSPs, então implementado Ce obtive os resultados matlabusando a função mex. Para calcular os coeficientes, use este comando em 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 )

DX
fonte