Algoritmo Goertzel - Goertzel algorithm

O algoritmo de Goertzel é uma técnica em processamento digital de sinais (DSP) para avaliação eficiente dos termos individuais da transformada discreta de Fourier (DFT). É útil em certas aplicações práticas, como o reconhecimento de tons de sinalização multifrequência de tom dual (DTMF) produzidos pelos botões do teclado de um telefone analógico tradicional . O algoritmo foi descrito pela primeira vez por Gerald Goertzel em 1958.

Como o DFT, o algoritmo Goertzel analisa um componente de frequência selecionável de um sinal discreto . Ao contrário dos cálculos DFT diretos, o algoritmo Goertzel aplica um único coeficiente de valor real a cada iteração, usando aritmética de valor real para sequências de entrada de valor real. Para cobrir um espectro completo, o algoritmo de Goertzel tem uma ordem de complexidade maior do que os algoritmos de transformada rápida de Fourier (FFT), mas para calcular um pequeno número de componentes de frequência selecionados, é mais eficiente numericamente. A estrutura simples do algoritmo Goertzel o torna adequado para pequenos processadores e aplicativos incorporados.

O algoritmo Goertzel também pode ser usado "ao contrário" como uma função de síntese sinusóide, que requer apenas 1 multiplicação e 1 subtração por amostra gerada.

O algoritmo

O cálculo principal do algoritmo de Goertzel tem a forma de um filtro digital e, por esta razão, o algoritmo é freqüentemente chamado de filtro de Goertzel . O filtro opera em uma seqüência de entrada em uma cascata de dois estágios com um parâmetro , dando a freqüência a ser analisada, normalizada em radianos por amostra. ${\ displaystyle x [n]}$${\ displaystyle \ omega _ {0}}$

O primeiro estágio calcula uma sequência intermediária : ${\ displaystyle s [n]}$

${\ displaystyle s [n] = x [n] +2 \ cos (\ omega _ {0}) s [n-1] -s [n-2].}$

(1)

O segundo estágio aplica o seguinte filtro a , produzindo a sequência de saída : ${\ displaystyle s [n]}$${\ displaystyle y [n]}$

${\ displaystyle y [n] = s [n] -e ^ {- j \ omega _ {0}} s [n-1].}$

(2)

O primeiro estágio de filtro pode ser observado como um filtro IIR de segunda ordem com uma estrutura de forma direta . Essa estrutura particular tem a propriedade de que suas variáveis ​​de estado internas sejam iguais aos valores de saída anteriores desse estágio. Os valores de entrada para são presumidos todos iguais a 0. Para estabelecer o estado do filtro inicial para que a avaliação possa começar na amostra , os estados do filtro são atribuídos a valores iniciais . Para evitar riscos de aliasing , a frequência é freqüentemente restrita ao intervalo de 0 a π (consulte o teorema de amostragem de Nyquist-Shannon ); usar um valor fora desse intervalo não é sem sentido, mas é equivalente a usar uma frequência com alias dentro desse intervalo, uma vez que a função exponencial é periódica com um período de 2π pol . ${\ displaystyle x [n]}$${\ displaystyle n <0}$${\ displaystyle x [0]}$${\ displaystyle s [-2] = s [-1] = 0}$${\ displaystyle \ omega _ {0}}$${\ displaystyle \ omega _ {0}}$

O filtro de segundo estágio pode ser observado como um filtro FIR , uma vez que seus cálculos não utilizam nenhuma de suas saídas anteriores.

Métodos de transformada Z podem ser aplicados para estudar as propriedades da cascata de filtros. A transformada Z do primeiro estágio do filtro dada na equação (1) é

${\ displaystyle {\ begin {align} {\ frac {S (z)} {X (z)}} & = {\ frac {1} {1-2 \ cos (\ omega _ {0}) z ^ { -1} + z ^ {- 2}}} \\ & = {\ frac {1} {(1-e ^ {+ j \ omega _ {0}} z ^ {- 1}) (1-e ^ {-j \ omega _ {0}} z ^ {- 1})}}. \ end {alinhado}}}$

(3)

A transformada Z do segundo estágio do filtro dada na equação (2) é

${\ displaystyle {\ frac {Y (z)} {S (z)}} = 1-e ^ {- j \ omega _ {0}} z ^ {- 1}.}$

(4)

A função de transferência combinada da cascata dos dois estágios do filtro é então

${\ displaystyle {\ begin {alinhado} {\ frac {S (z)} {X (z)}} {\ frac {Y (z)} {S (z)}} = {\ frac {Y (z) } {X (z)}} & = {\ frac {(1-e ^ {- j \ omega _ {0}} z ^ {- 1})} {(1-e ^ {+ j \ omega _ { 0}} z ^ {- 1}) (1-e ^ {- j \ omega _ {0}} z ^ {- 1})}} \\ & = {\ frac {1} {1-e ^ { + j \ omega _ {0}} z ^ {- 1}}}. \ end {alinhado}}}$

(5)

Isso pode ser transformado de volta em uma sequência equivalente no domínio do tempo, e os termos desenrolados de volta para o primeiro termo de entrada no índice : ${\ displaystyle n = 0}$

${\ displaystyle {\ begin {alinhados} y [n] & = x [n] + e ^ {+ j \ omega _ {0}} y [n-1] \\ & = \ sum _ {k = - \ infty} ^ {n} x [k] e ^ {+ j \ omega _ {0} (nk)} \\ & = e ^ {j \ omega _ {0} n} \ sum _ {k = 0} ^ {n} x [k] e ^ {- j \ omega _ {0} k} \ qquad {\ text {desde}} \ forall k <0, x [k] = 0. \ end {alinhado}}}$

(6)

Estabilidade numérica

Pode-se observar que os pólos da transformada Z do filtro estão localizados em e , em um círculo de raio unitário centrado na origem do plano da transformada Z complexa. Esta propriedade indica que o processo de filtro é marginalmente estável e vulnerável ao acúmulo de erros numéricos quando calculado usando aritmética de baixa precisão e sequências de entrada longas. Uma versão numericamente estável foi proposta por Christian Reinsch. ${\ displaystyle e ^ {+ j \ omega _ {0}}}$${\ displaystyle e ^ {- j \ omega _ {0}}}$

Cálculos DFT

Para o caso importante de cálculo de um termo DFT, as seguintes restrições especiais são aplicadas.

• A filtragem termina no índice , onde é o número de termos na sequência de entrada do DFT.${\ displaystyle n = N}$${\ displaystyle N}$
• As frequências escolhidas para a análise de Goertzel são restritas à forma especial

${\ displaystyle \ omega _ {0} = 2 \ pi {\ frac {k} {N}}.}$

(7)

• O número do índice que indica o "bin de frequência" do DFT é selecionado a partir do conjunto de números do índice${\ displaystyle k}$

${\ displaystyle k \ in \ {0,1,2, ..., N-1 \}.}$

(8)

Fazendo essas substituições na equação (6) e observando que o termo , equação (6) assume a seguinte forma: ${\ displaystyle e ^ {+ j2 \ pi k} = 1}$

${\ displaystyle y [N] = \ sum _ {n = 0} ^ {N} x [n] e ^ {- j2 \ pi {\ frac {nk} {N}}}.}$

(9)

Podemos observar que o lado direito da equação (9) é extremamente semelhante à fórmula de definição para o termo DFT , o termo DFT para o número índice , mas não exatamente o mesmo. O somatório mostrado na equação (9) requer termos de entrada, mas apenas os termos de entrada estão disponíveis ao avaliar um DFT. Um expediente simples, mas deselegante, é estender a sequência de entrada com mais um valor artificial . Podemos ver na equação (9) que o efeito matemático no resultado final é o mesmo que remover o termo da soma, entregando assim o valor DFT pretendido. ${\ displaystyle X [k]}$${\ displaystyle k}$${\ displaystyle N + 1}$${\ displaystyle N}$${\ displaystyle x [n]}$${\ displaystyle x [N] = 0}$${\ displaystyle x [N]}$

No entanto, existe uma abordagem mais elegante que evita a passagem extra do filtro. Da equação (1), podemos notar que quando o termo de entrada estendido é usado na etapa final, ${\ displaystyle x [N] = 0}$

${\ displaystyle s [N] = 2 \ cos (\ omega _ {0}) s [N-1] -s [N-2].}$

(10)

Assim, o algoritmo pode ser concluído da seguinte forma:

• encerrar o filtro IIR após o processamento do termo de entrada ,${\ displaystyle x [N-1]}$
• aplique a equação (10) para construir a partir dos resultados anteriores e ,${\ displaystyle s [N]}$${\ displaystyle s [N-1]}$${\ displaystyle s [N-2]}$
• aplique a equação (2) com o valor calculado e com o produzido pelo cálculo direto final do filtro.${\ displaystyle s [N]}$${\ displaystyle s [N-1]}$

As duas últimas operações matemáticas são simplificadas combinando-as algebricamente:

${\ displaystyle {\ begin {alinhados} y [N] & = s [N] -e ^ {- j2 \ pi {\ frac {k} {N}}} s [N-1] \\ & = (2 \ cos (\ omega _ {0}) s [N-1] -s [N-2]) - e ^ {- j2 \ pi {\ frac {k} {N}}} s [N-1] \ \ & = e ^ {j2 \ pi {\ frac {k} {N}}} s [N-1] -s [N-2]. \ end {alinhado}}}$

(11)

Observe que interromper as atualizações do filtro no final e aplicar imediatamente a equação (2) em vez da equação (11) perde as atualizações finais do estado do filtro, gerando um resultado com fase incorreta. ${\ displaystyle N-1}$

A estrutura de filtragem particular escolhida para o algoritmo Goertzel é a chave para seus cálculos DFT eficientes. Podemos observar que apenas um valor de saída é usado para calcular o DFT, portanto, os cálculos para todos os outros termos de saída são omitidos. Como o filtro FIR não é calculado, os cálculos do estágio IIR , etc. podem ser descartados imediatamente após a atualização do estado interno do primeiro estágio. ${\ displaystyle y [N]}$${\ displaystyle s [0], s [1]}$

Isso parece deixar um paradoxo: para completar o algoritmo, o estágio do filtro FIR deve ser avaliado uma vez usando as duas saídas finais do estágio do filtro IIR, enquanto para eficiência computacional a iteração do filtro IIR descarta seus valores de saída. É aqui que as propriedades da estrutura de filtro de forma direta são aplicadas. As duas variáveis ​​de estado internas do filtro IIR fornecem os dois últimos valores da saída do filtro IIR, que são os termos necessários para avaliar o estágio do filtro FIR.

Formulários

Termos do espectro de potência

Examinando a equação (6), uma passagem final do filtro IIR para calcular o termo usando um valor de entrada suplementar aplica um multiplicador complexo de magnitude 1 ao termo anterior . Conseqüentemente, e representam potência de sinal equivalente. É igualmente válido aplicar a equação (11) e calcular a potência do sinal a partir do termo ou aplicar a equação (2) e calcular a potência do sinal a partir do termo . Ambos os casos levam à seguinte expressão para a potência do sinal representada pelo termo DFT : ${\ displaystyle y [N]}$${\ displaystyle x [N] = 0}$${\ displaystyle y [N-1]}$${\ displaystyle y [N]}$${\ displaystyle y [N-1]}$${\ displaystyle y [N]}$${\ displaystyle y [N-1]}$${\ displaystyle X [k]}$

${\ displaystyle {\ begin {alinhado} X [k] \, X '[k] & = y [N] \, y' [N] = y [N-1] \, y '[N-1] \ \ & = s ^ {2} [N-1] + s ^ {2} [N-2] -2 \ cos \ left (2 \ pi {\ frac {k} {N}} \ right) \, s [N-1] \, s [N-2]. \ End {alinhado}}}$

(12)

No pseudocódigo abaixo, as variáveis spreve sprev2armazenam temporariamente o histórico de saída do filtro IIR, enquanto x[n]é um elemento indexado da matriz x , que armazena a entrada.

Nterms defined here
Kterm selected here
ω = 2 × π × Kterm / Nterms;
cr := cos(ω)
ci := sin(ω)
coeff := 2 × cr

sprev := 0
sprev2 := 0
for each index n in range 0 to Nterms-1 do
    s := x[n] + coeff × sprev - sprev2
    sprev2 := sprev
    sprev := s
end

power := sprev2 × sprev2 + sprev × sprev - coeff × sprev × sprev2

É possível organizar os cálculos de forma que as amostras de entrada sejam entregues individualmente a um objeto de software que mantém o estado do filtro entre as atualizações, com o resultado de potência final acessado após o outro processamento ser feito.

Termo DFT único com aritmética de valor real

O caso de dados de entrada com valor real surge frequentemente, especialmente em sistemas embarcados onde os fluxos de entrada resultam de medições diretas de processos físicos. Comparando com a ilustração na secção anterior, quando os dados de entrada são reais de valor, as variáveis de estado internos de filtro spreve sprev2pode ser observado também ser de valor real, consequentemente, não aritmética complexo é necessária na primeira etapa IIR. A otimização para aritmética de valor real normalmente é tão simples quanto aplicar os tipos de dados de valor real apropriados para as variáveis.

Após os cálculos usando o termo de entrada e as iterações do filtro serem finalizados, a equação (11) deve ser aplicada para avaliar o termo DFT. O cálculo final usa aritmética de valor complexo, mas isso pode ser convertido em aritmética de valor real separando termos reais e imaginários: ${\ displaystyle x [N-1]}$

${\ displaystyle {\ begin {alinhados} c_ {r} & = \ cos (2 \ pi {\ tfrac {k} {N}}), \\ c_ {i} & = \ sin (2 \ pi {\ tfrac {k} {N}}), \\ y [N] & = c_ {r} s [N-1] -s [N-2] + jc_ {i} s [N-1]. \ end {alinhado }}}$

(13)

Comparando com a aplicação do espectro de potência, a única diferença é o cálculo usado para terminar:

(Same IIR filter calculations as in the signal power implementation)
XKreal = sprev * cr - sprev2;
XKimag = sprev * ci;

Detecção de fase

Esta aplicação requer a mesma avaliação do termo DFT , conforme discutido na seção anterior, usando um fluxo de entrada de valor real ou complexo. Então, a fase do sinal pode ser avaliada como ${\ displaystyle X [k]}$

${\ displaystyle \ phi = \ tan ^ {- 1} {\ frac {\ Im (X [k])} {\ Re (X [k])}},}$

(14)

tomando as precauções adequadas para singularidades, quadrantes e assim por diante ao calcular a função tangente inversa.

Sinais complexos em aritmética real

Uma vez que sinais complexos se decompõem linearmente em partes reais e imaginárias, o algoritmo Goertzel pode ser calculado em aritmética real separadamente sobre a sequência de partes reais, cedendo , e sobre a sequência de partes imaginárias, cedendo . Depois disso, os dois resultados parciais de valor complexo podem ser recombinados: ${\ displaystyle y _ {\ text {r}} [n]}$${\ displaystyle y _ {\ text {i}} [n]}$

${\ displaystyle y [n] = y _ {\ text {r}} [n] + jy _ {\ text {i}} [n].}$

(15)

Complexidade computacional

• De acordo com a teoria da complexidade computacional , o cálculo de um conjunto de termos DFT usando aplicativos do algoritmo Goertzel em um conjunto de dados com valores com um "custo por operação" de tem complexidade .${\ displaystyle M}$${\ displaystyle M}$${\ displaystyle N}$${\ displaystyle K}$ ${\ displaystyle O (KNM)}$

Para calcular um único bin DFT para uma sequência de entrada complexa de comprimento , o algoritmo Goertzel requer multiplicações e adições / subtrações dentro do loop, bem como 4 multiplicações e 4 adições / subtrações finais, para um total de multiplicações e adições / subtrações. Isso é repetido para cada uma das frequências.${\ displaystyle X (f)}$${\ displaystyle N}$${\ displaystyle 2N}$${\ displaystyle 4 \ N}$${\ displaystyle 2N + 4}$${\ displaystyle 4N + 4}$${\ displaystyle M}$

• Em contraste, usar um FFT em um conjunto de dados com valores tem complexidade .${\ displaystyle N}$${\ displaystyle O (KN \ log N)}$

Isso é mais difícil de aplicar diretamente porque depende do algoritmo FFT usado, mas um exemplo típico é um FFT radical-2, que requer multiplicações e adições / subtrações por bin DFT , para cada um dos bins.${\ displaystyle 2 \ log _ {2} (N)}$${\ displaystyle 3 \ log _ {2} (N)}$${\ displaystyle N}$

Nas expressões de ordem de complexidade, quando o número de termos calculados é menor que , a vantagem do algoritmo de Goertzel é clara. Mas como o código FFT é comparativamente complexo, o fator de "custo por unidade de trabalho" costuma ser maior para um FFT, e a vantagem prática favorece o algoritmo de Goertzel até mesmo para várias vezes maior do que . ${\ displaystyle M}$${\ displaystyle \ log N}$${\ displaystyle K}$${\ displaystyle M}$${\ displaystyle \ log _ {2} (N)}$

Como regra geral para determinar se um FFT radix-2 ou um algoritmo de Goertzel é mais eficiente, ajuste o número de termos no conjunto de dados para cima para a potência exata mais próxima de 2, chamando isso , e o algoritmo de Goertzel é provável ser mais rápido se ${\ displaystyle N}$${\ displaystyle N_ {2}}$

${\ displaystyle M \ leq {\ frac {5N_ {2}} {6N}} \ log _ {2} (N_ {2})}$

As implementações de FFT e as plataformas de processamento têm um impacto significativo no desempenho relativo. Algumas implementações de FFT executam cálculos de números complexos internos para gerar coeficientes em tempo real, aumentando significativamente seu "custo K por unidade de trabalho". Os algoritmos FFT e DFT podem usar tabelas de valores de coeficientes pré-calculados para melhor eficiência numérica, mas isso requer mais acessos aos valores de coeficientes armazenados em buffer na memória externa, o que pode levar a uma maior contenção de cache que contraria algumas das vantagens numéricas.

Ambos os algoritmos ganham aproximadamente um fator de 2 de eficiência ao usar dados de entrada de valor real em vez de dados de valor complexo. No entanto, esses ganhos são naturais para o algoritmo Goertzel, mas não serão alcançados para o FFT sem o uso de certas variantes de algoritmo especializadas para transformar dados de valor real .

Veja também

• Algoritmo FFT de Bluestein (chirp-Z)
• Mudança de frequência (FSK)
• Phase-shift keying (PSK)

Referências

Leitura adicional

• Proakis, JG; Manolakis, DG (1996), Digital Signal Processing: Principles, Algorithms, and Applications , Upper Saddle River, NJ: Prentice Hall, pp. 480-481, Bibcode : 1996dspp.book ..... P

links externos

• Algoritmo Goertzel na Máquina Wayback (arquivado em 28/06/2018)
• Um algoritmo DSP para análise de frequência
• O Algoritmo Goertzel de Kevin Banks