Obtenir la valeur de crête d'un signal si sa fréquence se situe entre deux centres de bin

12

Veuillez supposer ce qui suit:

  • La fréquence du fondamental d'un signal a été estimée à l'aide de la FFT et de certaines méthodes d'estimation de fréquence et se situe entre deux centres de bin
  • La fréquence d'échantillonnage est fixe
  • L'effort de calcul n'est pas un problème

Connaissant la fréquence, quelle est la manière la plus précise d'estimer la valeur de crête correspondante des signaux fondamentaux?

Une façon pourrait être de mettre à zéro le signal temporel pour augmenter la résolution FFT de sorte que le centre du bac soit plus proche de la fréquence estimée. Dans ce scénario, un point dont je ne suis pas sûr est de savoir si je peux mettre à zéro autant que je veux ou s'il y a des inconvénients à le faire. Un autre est le centre du bac que je devrais sélectionner après un remplissage nul comme celui dont j'obtiens la valeur de crête (car on ne peut pas atteindre la fréquence d'intérêt exactement, même après un remplissage nul).

Cependant, je me demande également s'il existe une autre méthode qui peut donner de meilleurs résultats, par exemple un estimateur qui utilise les valeurs de crête des deux centres de cellules environnants pour estimer la valeur de crête à la fréquence d'intérêt.

lR8n6i
la source
2
zéro padding avant FFT est un moyen. Une autre consiste à appliquer une fonction de fenêtre adaptée à vos neads. La fenêtre supérieure plate a été conçue exactement dans ce but. Bien sûr, si vous connaissez déjà la fréquence exactement et que vous êtes intéressé par un seul amplutide, il existe probablement des moyens moins chers de le faire qu'une FFT.
sellibitze
1
pas de remplissage nul requis: une simple interpolation parabolique (avec 3 points: imax-1, imax, imax + 1, où imaxest le pic FFT) vous donnera des résultats précis
Basj
Assurez-vous que la fonction d'interpolation correspond à la fonction de fenêtre. Flat-top est trivial, sinon vous voulez une paire correspondante (par exemple fenêtre rectangulaire + interpolation sinc, fenêtre gaussienne + interpolation gaussienne, etc.)
finnw
@CedronDawg cette question et ses réponses sont liées (mais pas les mêmes) à votre formule de fréquence exacte. Peut-être que vous pouvez le trouver intéressant.
Fat32

Réponses:

5

Le premier algorithme qui me vient à l'esprit est l' algorithme de Goertzel . Cet algorithme suppose généralement que la fréquence d'intérêt est un multiple entier de la fréquence fondamentale. Cependant, cet article applique l'algorithme (généralisé) au cas qui vous intéresse.


Un autre problème est que le modèle de signal est incorrect. Il utilise 2*%pi*(1:siglen)*(Fc/siglen). Il devrait utiliser 2*%pi*(0:siglen-1)*(Fc/siglen)pour que la phase sorte correctement.

Je pense aussi qu'il y a un problème avec la fréquence Fc=21.3très basse. Les signaux à basse fréquence à valeur réelle ont tendance à présenter un biais lorsqu'il s'agit de problèmes d'estimation de phase / fréquence.

J'ai également essayé une recherche de grille grossière pour l'estimation de phase, et elle donne la même réponse que l'algorithme de Goertzel.

Vous trouverez ci-dessous un graphique qui montre le biais dans les deux estimations (Goertzel: bleu, grossier: rouge) pour deux fréquences différentes: Fc=21.3(solide) et Fc=210.3(en pointillés). Comme vous pouvez le voir, le biais pour la fréquence plus élevée est beaucoup moins.

x2π

entrez la description de l'image ici

Peter K.
la source
Je viens de tester le code de l'algorithme Goerzel basé sur le papier. En utilisant la valeur DTFT de sortie, le pic peut être obtenu très précisément. Cependant, il y a un facteur d'échelle exactement de 1000. Donc, si le pic d'origine est de 1 234, après Goerzel, ce sera 1234. Quelqu'un sait-il d'où cela pourrait provenir?
lR8n6i
A fait quelques recherches dans l'intervalle. Cela a probablement à voir avec la mise à l'échelle de l'amplitude: mise à l'échelle de l'amplitude du domaine temporel = coefficient du domaine fréquentiel * 2 / N, où N est la longueur du signal. Cette hypothèse est-elle juste?
lR8n6i
Salut! Je viens de découvrir qu'en utilisant l'algorithme de Goertzel, l'amplitude au coefficient complexe résultant est très précise, mais la phase est complètement fausse. Quelqu'un a-t-il une idée d'où cela pourrait provenir? Par "phase", j'entends le décalage de phase spécifié dans la fondamentale du signal d'origine.
lR8n6i
1
sin(ω0t+ϕ)j2[ejϕδ~(ω+ω0+2πk)e+jϕδ~(ωω0+2πk)]π/2
4

Si vous êtes prêt à utiliser plusieurs bacs FFT voisins, pas seulement 2, l'interpolation Sinc fenêtrée entre les résultats de bac complexes peut produire une estimation très précise, en fonction de la largeur de la fenêtre.

L'interpolation Windinc Sinc se trouve couramment dans les suréchantillonneurs audio de haute qualité, de sorte que les articles sur ce sujet auront des formules d'interpolation appropriées avec analyse des erreurs.

hotpaw2
la source
Merci pour le commentaire. Je vais également essayer cette approche.
lR8n6i
4

sin(πx)(πx)

[1] JL Flanagan et RM Golden, «Phase vocoder», Bell Systems Technical Journal, vol. 45, p. 1493-1509, 1966.

[2] K. Dressler, «Extraction sinusoïdale utilisant une implémentation ef fi cace d'une FFT multi-résolution», dans Proc. 9th Int. Conf. sur les effets audio numériques (DAFx-06), Montréal, Canada, sept. 2006, pp. 247–252.

ederwander
la source
Salut! Merci beaucoup pour tous vos commentaires. J'ai étendu mon code (voir ci-dessous) pour combiner le filtre de Goertzel avec une interpolation de pic parabolique pour obtenir la phase. Cependant, les résultats ne sont toujours pas précis (+ - 3-4deg). Est-ce aussi proche que possible ou y a-t-il des erreurs de compréhension ou de codage?
lR8n6i
3

J'ai eu beaucoup de difficultés avec ce problème il y a quelques années.

J'ai posté cette question:

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

J'ai fini par faire les calculs à partir de zéro et j'ai posté une réponse à ma propre question.

Je suis surpris de n'avoir pu trouver aucune exposition similaire sur Internet.

Je posterai à nouveau la réponse ici; notez que le code est conçu pour un scénario dans lequel je chevauche ma fenêtre FFT de 4x.

π


Ce puzzle prend deux clés pour le déverrouiller.

Graphique 3.3:

entrez la description de l'image ici

Graphique 3.4:

entrez la description de l'image ici

Code:

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
la source
Vous interpolez la fréquence, alors que l'OP connaît la fréquence et veut interpoler l'amplitude.
finnw
2

Ce code python vous donnera un résultat très précis (je l'ai utilisé pour beaucoup de notes de musique et obtenu des erreurs inférieures à 0,01% de demi-ton) avec interpolation parabolique (méthode utilisée avec succès par McAulay Quatieri, Serra, etc. en harmonique + résiduel techniques de séparation)

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
la source
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

Les fréquences que vous traitez (21,3 Hz échantillonnés à 8 kHz) sont très basses. Parce que ce sont des signaux à valeur réelle, ils présenteront un biais dans l'estimation de phase pour ** n'importe quelle ** fréquence.

Cette image montre un tracé du biais ( phase_est - phase_orig) pour Fc = 210.3;(en rouge) par rapport au biais pour Fc = 21.3;. Comme vous pouvez le voir, le décalage est beaucoup plus important pour le 21.3cas.

Une autre option consiste à réduire votre taux d'échantillonnage. La courbe verte montre le biais pour Fs = 800au lieu de 8000.

entrez la description de l'image ici

lR8n6i
la source
1
Merci pour la mise à jour! Voir mon intrigue; Je pense toujours que tout estimateur de phase aura un biais pour cette basse fréquence. Une façon de le contourner est d'utiliser la fréquence connue (si elle est connue!) Pour corriger le biais d'estimation de phase à travers une table de correspondance. Mais vous devrez faire attention: le biais changera avec la fréquence. Une autre façon de le faire sera de réduire votre taux d'échantillonnage.
Peter K.
1
Merci aussi! Cependant, si vous utilisez Fs = 8000 Hz et Fc = 210 au lieu de 210,3, le biais semble encore pire. Une idée d'où cela pourrait venir?
lR8n6i
1
Erk! Aucune idée. FWIW, l'estimateur Geortzel n'a pas de problèmes: goertzel = atan(imag(y(2)),real(y(2)))*180/%pi + 90;. :-) Va creuser un peu plus. Surveillez cet endroit.
Peter K.
1
L'interpolation parabolique ne fait pas ce que vous pensez qu'elle fait. En particulier, si vous remplacez votre calcul ppar 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;alors vous obtenez BEAUCOUP de meilleures réponses --- même pour Fc=210. Je ne suis pas du tout sûr que la version actuelle de pvous donnera quelque chose de sensé. La formule d'interpolation est pour l'interpolation de l'AMPLITUDE d'une parabole, mais pinterpole la phase qui est juste ... bizarre.
Peter K.
1
Tout cela est OK, SAUF que l'emplacement du pic ( p = (yp1 - ym1)/(2*(2*y0 - yp1 - ym1))) sera incorrect parfois si vous utilisez des PHASES au lieu des amplitudes. En effet, les phases peuvent sauter autour de la limite de +/- 180 degrés. Tout ce qui est nécessaire pour le corriger pour la phase est de changer cette ligne pour mon p2calcul ci-dessus.
Peter K.