Obtenha o valor de pico de um sinal se sua frequência estiver entre dois centros de compartimento

12

Suponha o seguinte:

  • A frequência fundamental de um sinal foi estimada usando FFT e alguns métodos de estimativa de frequência e fica entre dois centros de bin
  • A frequência de amostragem é fixa
  • Esforço computacional não é um problema

Conhecendo a frequência, qual é a maneira mais precisa de estimar o valor de pico correspondente dos sinais fundamentais?

Uma maneira pode ser zerar o sinal de tempo para aumentar a resolução da FFT, de modo que o centro da bandeja esteja mais próximo da frequência estimada. Nesse cenário, um ponto sobre o qual não tenho certeza é se consigo zerar tanto quanto quero ou se há algumas desvantagens em fazê-lo. Outro é o centro de posição em que eu devo selecionar após o preenchimento zero como o de onde estou obtendo o valor de pico (porque não é possível atingir exatamente a frequência de interesse, mesmo após o preenchimento zero).

No entanto, também estou me perguntando se existe outro método que possa fornecer melhores resultados, por exemplo, um estimador que use os valores de pico dos dois centros de caixa ao redor para estimar o valor de pico na frequência de interesse.

lR8n6i
fonte
2
o preenchimento zero antes da FFT é uma maneira. Outra é aplicar uma função de janela adequada para seus neads. A janela superior plana foi projetada exatamente para esse fim. Obviamente, se você já conhece exatamente a frequência e está interessado em apenas um amplutídeo, provavelmente existem maneiras mais baratas de fazer isso do que uma FFT.
sellibitze
1
não é necessário preenchimento zero: interpolação parabólica simples (com 3 pontos: imax-1, imax, imax + 1, onde imaxestá o pico da FFT) fornecerá resultados precisos
Basj
Verifique se a função de interpolação corresponde à função da janela. Flat-top é trivial, caso contrário, você quer um par de correspondência (por exemplo janela + interpolação sinc retangular, janela Gaussian + gaussian interpolação etc.)
finnw
@CedronDawg esta pergunta e suas respostas estão relacionadas (mas não são as mesmas) com a sua fórmula exata de frequência. Pode ser que você possa achar interessante.
precisa

Respostas:

5

O primeiro algoritmo que vem à mente é o algoritmo de Goertzel . Esse algoritmo geralmente assume que a frequência de interesse é um múltiplo inteiro da frequência fundamental. No entanto, este artigo aplica o algoritmo (generalizado) ao caso em que você está interessado.


Outro problema é que o modelo de sinal está incorreto. Ele usa 2*%pi*(1:siglen)*(Fc/siglen). Deve ser usado 2*%pi*(0:siglen-1)*(Fc/siglen)para que a fase saia corretamente.

Eu também acho que há um problema com a frequência Fc=21.3muito baixa. Sinais com valor real de baixa frequência tendem a exibir viés quando se trata de problemas de estimativa de fase / frequência.

Eu também tentei uma pesquisa de grade grossa para a estimativa de fase, e ela fornece a mesma resposta que o algoritmo de Goertzel.

Abaixo está um gráfico que mostra o viés nas duas estimativas (Goertzel: azul, Grosso: vermelho) para duas frequências diferentes: Fc=21.3(sólida) e Fc=210.3(tracejada). Como você pode ver, o viés para a frequência mais alta é muito menor.

O gráfico eixo é a fase inicial que muda de 0 para 2 π .x2π

insira a descrição da imagem aqui

Peter K.
fonte
Acabei de testar o código do algoritmo de Goerzel com base no artigo. Usando o valor DTFT de saída, o pico pode ser obtido com muita precisão. No entanto, existe um fator de escala de exatamente 1000. Portanto, se o pico original for 1.234, depois de Goerzel será 1234. Alguém sabe de onde isso pode vir?
LR8n6i
Fiz algumas pesquisas nesse meio tempo. Provavelmente, isso tem a ver com a escala de amplitude: amplitude no domínio do tempo de escala = coeficiente no domínio da frequência * 2 / N, em que N é o comprimento do sinal. Essa suposição está certa?
LR8n6i
Oi! Acabei de descobrir que, usando o algoritmo de Goertzel, a amplitude no coeficiente complexo resultante é muito precisa, mas a fase está completamente errada. Alguém tem uma idéia de onde isso pode vir? Por "fase", quero dizer o atraso de fase especificado no fundamental do sinal original.
LR8n6i
1
sin(ω0t+ϕ)j2[ejϕδ~(ω+ω0+2πk)e+jϕδ~(ωω0+2πk)]π/2
4

Se você estiver disposto a usar vários compartimentos FFT vizinhos, não apenas 2, a interpolação Sinc com janela entre os resultados complexos da bandeja pode produzir uma estimativa muito precisa, dependendo da largura da janela.

A interpolação Windinc Sinc é comumente encontrada em amplificadores de áudio de alta qualidade, portanto, os trabalhos sobre esse assunto terão fórmulas de interpolação adequadas com a análise de erros.

hotpaw2
fonte
Obrigado pelo comentário. Vou tentar esta abordagem também.
LR8n6i
4

pecado(πx)(πx)

[1] JL Flanagan e RM Golden, "Phase vocoder", Bell Systems Technical Journal, vol. 45, pp. 1493–1509, 1966.

[2] K. Dressler, “Extração sinusoidal usando uma implementação eficiente de uma FFT de alta resolução”, no Proc. 9ª Int. Conf. sobre efeitos de áudio digital (DAFx-06), Montreal, Canadá, setembro de 2006, pp. 247–252.

Ederwander
fonte
Oi! Muito obrigado por todos os seus comentários. Estendi meu código (veja abaixo) para combinar o filtro Goertzel com a interpolação de pico parabólico para obter a fase. No entanto, os resultados ainda não são precisos (+ - 3-4 graus). É o mais próximo possível ou há erros no entendimento ou na codificação?
LR8n6i
3

Eu tive muita dificuldade com esse problema exato alguns anos atrás.

Eu postei esta pergunta:

/programming/4633203/extracting-precise-frequencies-from-fft-bins-using-phase-change-between-frames

Acabei fazendo os cálculos do zero e postei uma resposta para minha própria pergunta.

Estou surpreso por não ter encontrado nenhuma exposição semelhante na Internet.

Vou postar a resposta novamente aqui; observe que o código foi projetado para um cenário em que estou sobrepondo minha janela FFT em 4x.

π


Este quebra-cabeça leva duas chaves para desbloqueá-lo.

Gráfico 3.3:

insira a descrição da imagem aqui

Gráfico 3.4:

insira a descrição da imagem aqui

Código:

for (int k = 0; k <= fftFrameSize/2; k++) 
{
    // compute magnitude and phase 
    bins[k].mag = 2.*sqrt(fftBins[k].real*fftBins[k].real + fftBins[k].imag*fftBins[k].imag);
    bins[k].phase = atan2(fftBins[k].imag, fftBins[k].real);

    // Compute phase difference Δϕ fo bin[k]
    double deltaPhase;
    {
        double measuredPhaseDiff = bins[k].phase - gLastPhase[k];
        gLastPhase[k] = bins[k].phase;

        // Subtract expected phase difference <-- FIRST KEY
        // Think of a single wave in a 1024 float frame, with osamp = 4
        //   if the first sample catches it at phase = 0, the next will 
        //   catch it at pi/2 ie 1/4 * 2pi
        double binPhaseExpectedDiscrepancy = M_TWOPI * (double)k / (double)osamp;
        deltaPhase = measuredPhaseDiff - binPhaseExpectedDiscrepancy;

        // Wrap delta phase into [-Pi, Pi) interval 
        deltaPhase -= M_TWOPI * floor(deltaPhase / M_TWOPI + .5);
    }

    // say sampleRate = 40K samps/sec, fftFrameSize = 1024 samps in FFT giving bin[0] thru bin[512]
    // then bin[1] holds one whole wave in the frame, ie 44 waves in 1s ie 44Hz ie sampleRate / fftFrameSize
    double bin1Freq = (double)sampleRate / (double)fftFrameSize;
    bins[k].idealFreq = (double)k * bin1Freq;

    // Consider Δϕ for bin[k] between hops.
    // write as 2π / m.
    // so after m hops, Δϕ = 2π, ie 1 extra cycle has occurred   <-- SECOND KEY
    double m = M_TWOPI / deltaPhase;

    // so, m hops should have bin[k].idealFreq * t_mHops cycles.  plus this extra 1.
    // 
    // bin[k].idealFreq * t_mHops + 1 cycles in t_mHops seconds 
    //   => bins[k].actualFreq = bin[k].idealFreq + 1 / t_mHops
    double tFrame = fftFrameSize / sampleRate;
    double tHop = tFrame / osamp;
    double t_mHops = m * tHop;

    bins[k].freq = bins[k].idealFreq + 1. / t_mHops;
}
P i
fonte
Você está interpolando a frequência, enquanto o OP conhece a frequência e deseja interpolar a amplitude.
finnw
2

Este código python fornecerá um resultado muito preciso (usei-o para muitas notas musicais e obtive erros inferiores a 0,01% do semitom) com interpolação parabólica (método usado com sucesso por McAulay Quatieri, Serra, etc. em harmônico + residual técnicas de separação)

import matplotlib.pyplot as plt
import numpy as np
from scipy.io.wavfile import read
from scipy.fftpack import fft, ifft
import math

(fs, x) = read('test.wav')
if (len(x.shape) == 2):    # if stereo we keep left channel only
 x = x[:,1]

n=x.size
freq = np.arange(n)*1.0/n*fs 
xfft = abs(fft(x))

imax=np.argmax(xfft)  
p=1.0/2*(xfft[imax-1]/xfft[imax]-xfft[imax+1]/xfft[imax])/(xfft[imax-1]/xfft[imax]-2+xfft[imax+1]/xfft[imax])   # parabolic interpolation 
print 'Frequence detectee avec interpolation parabolique :',(imax+p)*1.0/n*fs, 'Hz'
Basj
fonte
1
clear all
clc

for phase_orig = 0:pi/18:pi,

%% Specify and generate signal
Amp = 1;                     % Amplitude of signal
Fs = 8000;                   % samples per second
dt = 1/Fs;                   % seconds per sample
Fc = 21.3;                   % Hz
StopTime = 0.25;             % seconds
t = (0:dt:StopTime-dt)';     % seconds

siglen = length(t);
sig = Amp * 1.5 * sin(2*pi*(0:siglen-1)*(Fc/siglen) + phase_orig) + 1.5 * Amp * sin(2*pi*(0:siglen-1)*(Fc/siglen) * 3) ...
  + 1.5 * Amp * sin(2*pi*(0:siglen-1)*(Fc/siglen) * 5)+ 0.3 * Amp * sin(2*pi*(0:siglen-1)*(Fc/siglen) * 7) ...
  + 1.3 * Amp * sin(2*pi*(0:siglen-1)*(Fc/siglen) * 9)+ 1.4 * Amp * sin(2*pi*(0:siglen-1)*(Fc/siglen) * 11);

%% Estimate the peak value of the signals fundamental using Goertzel algorithm
peak = 0;
indvec = [Fc-1 Fc Fc+1];

% Check the input data
if ~isvector(sig) || isempty(sig)
  error('X must be a nonempty vector')
end

if ~isvector(indvec) || isempty(indvec)
  error('INDVEC must be a nonempty vector')
end
if ~isreal(indvec)
  error('INDVEC must contain real numbers')
end

% forcing x to be column
sig = reshape(sig,siglen,1);

% initialization
no_freq = length(indvec); %number of frequencies to compute
y = zeros(no_freq,1); %memory allocation for the output coefficients

% Computation via second-order system
% loop over the particular frequencies
for cnt_freq = 1:no_freq
  %for a single frequency:
  %a/ precompute the constants
  pik_term = 2*pi*(indvec(cnt_freq))/(siglen);
  cos_pik_term2 = cos(pik_term) * 2;
  cc = exp(-1i*pik_term); % complex constant
  %b/ state variables
  s0 = 0;
  s1 = 0;
  s2 = 0;
  %c/ 'main' loop
  for ind = 1:siglen-1 %number of iterations is (by one) less than the length of signal
    %new state
    s0 = sig(ind) + cos_pik_term2 * s1 - s2;  % (*)
    %shifting the state variables
    s2 = s1;
    s1 = s0;
  end
  %d/ final computations
  s0 = sig(siglen) + cos_pik_term2 * s1 - s2; %correspond to one extra performing of (*)
  y(cnt_freq) = s0 - s1*cc; %resultant complex coefficient

  %complex multiplication substituting the last iterationA
  %and correcting the phase for (potentially) non-integer valued
  %frequencies at the same time
  y(cnt_freq) = y(cnt_freq) * exp(-1i*pik_term*(siglen-1));
end

  % perfom amplitude scaling
  peak = abs(y(2)) * 2 / siglen

% perform parabolic interpolation to get the phase estimate
phase_orig=phase_orig*180/pi
ym1 = angle(unwrap(y(1)));
y0 = angle(unwrap(y(2)));
yp1 = angle(unwrap(y(3)));

p = (yp1 - ym1)/(2*(2*y0 - yp1 - ym1)); 
phase = y0 - 0.25*(ym1-yp1)*p;
phase_est = phase * 180/pi + 90;
phase_est = mod(phase_est+180,360)-180
end

As frequências com as quais você está lidando (21,3Hz amostrados em 8kHz) são muito baixas. Por serem sinais de valor real, exibirão um viés na estimativa de fase para ** qualquer ** frequência.

Esta imagem mostra um gráfico do viés ( phase_est - phase_orig) para Fc = 210.3;(em vermelho) versus o viés para Fc = 21.3;. Como você pode ver, o deslocamento é muito mais significativo para o21.3 caso.

Outra opção é reduzir sua taxa de amostragem. A curva verde mostra o viés paraFs = 800 vez de 8000.

insira a descrição da imagem aqui

lR8n6i
fonte
1
Obrigado pela atualização! Veja meu enredo; Eu ainda acho que qualquer estimador de fase terá um viés para essa frequência baixa. Uma maneira de contornar isso é usar a frequência conhecida (se for conhecida!) Para corrigir o viés da estimativa de fase através de uma tabela de consulta. Mas você precisará ter cuidado: o viés mudará com a frequência. Outra maneira de fazer isso será reduzir sua taxa de amostragem.
Peter K.
1
Obrigada também! No entanto, se você estiver usando Fs = 8000 Hz e Fc = 210 em vez de 210,3, o viés parecerá ainda pior. Alguma idéia de onde isso poderia vir?
LR8n6i
1
Erk! Nenhuma idéia. FWIW, o estimador Geortzel não tem problemas: goertzel = atan(imag(y(2)),real(y(2)))*180/%pi + 90;. :-) Vai cavar um pouco mais. Assista esse espaço.
Peter K.
1
A interpolação parabólica não está fazendo o que você pensa que está fazendo. Em particular, se você substituir seu cálculo ppor p2 = (abs(y(3)) - abs(y(1)))/(2*(2*abs(y(2)) - abs(y(3)) - abs(y(1)))); phase2 = y0 - 0.25*(ym1-yp1)*p2;, obterá MUITAS melhores respostas - mesmo para Fc=210. Não tenho certeza de que a versão atual do plhe dará algo sensato. A fórmula de interpolação é para interpolação da AMPLITUDE de uma parábola, mas pestá interpolando a fase que é ... estranha.
Peter K.
1
Tudo isso está p = (yp1 - ym1)/(2*(2*y0 - yp1 - ym1))correto , EXCETO que a localização do pico ( ) estará incorreta algumas vezes se você estiver usando PHASES em vez de amplitudes. Isso ocorre porque as fases podem pular em torno do limite de +/- 180 graus. Tudo o que é necessário para corrigi-lo para a fase é alterar essa linha para o meu p2cálculo acima.
Peter K.