Ermitteln Sie den Spitzenwert eines Signals, wenn seine Frequenz zwischen zwei Bin-Zentren liegt

12

Bitte nehmen Sie folgendes an:

  • Die Frequenz der Grundschwingung eines Signals wurde unter Verwendung der FFT und einiger Frequenzschätzungsverfahren geschätzt und liegt zwischen zwei Bin-Zentren
  • Die Abtastfrequenz ist fest eingestellt
  • Rechenaufwand ist kein Thema

Wie lässt sich bei Kenntnis der Frequenz der entsprechende Spitzenwert der Grundsignale am genauesten abschätzen?

Eine Möglichkeit könnte darin bestehen, das Zeitsignal auf Null zu setzen, um die FFT-Auflösung so zu erhöhen, dass die Bin-Mitte näher an der geschätzten Frequenz liegt. In diesem Szenario bin ich mir nicht sicher, ob ich so viel Zero-Pad ausführen kann, wie ich möchte, oder ob dies einige Nachteile hat. Ein anderes ist, welches Bin-Zentrum ich nach dem Null-Auffüllen als das auswählen soll, von dem ich den Spitzenwert erhalte (weil man die interessierende Frequenz möglicherweise auch nach dem Null-Auffüllen nicht genau trifft).

Ich frage mich jedoch auch, ob es eine andere Methode gibt, die möglicherweise bessere Ergebnisse liefert, z. B. einen Schätzer, der die Spitzenwerte der beiden umgebenden Bin-Zentren verwendet, um den Spitzenwert bei der interessierenden Frequenz zu schätzen.

IR8n6i
quelle
2
Null-Auffüllung vor FFT ist eine Möglichkeit. Eine andere Möglichkeit besteht darin, eine Fensterfunktion anzuwenden, die für Ihre Neads geeignet ist. Das flache obere Fenster wurde genau für diesen Zweck entworfen. Wenn Sie die Frequenz bereits genau kennen und nur an einem Amplitudenverstärker interessiert sind, gibt es wahrscheinlich günstigere Möglichkeiten als eine FFT.
Sellibitze
1
Keine Null-Polsterung erforderlich: Einfache parabolische Interpolation (mit 3 Punkten: imax-1, imax, imax + 1, wo imaxist der FFT-Peak) gibt Ihnen genaue Ergebnisse
Basis
Stellen Sie sicher, dass die Interpolationsfunktion mit der Fensterfunktion übereinstimmt. Flat-Top ist trivial, ansonsten möchten Sie ein passendes Paar (z. B. rechteckiges Fenster + Sinc-Interpolation, Gauß-Fenster + Gauß-Interpolation usw.)
finnw
@CedronDawg Diese Frage und ihre Antworten sind mit Ihrer genauen Frequenzformel verknüpft (aber nicht identisch). Vielleicht finden Sie es interessant.
Fat32

Antworten:

5

Der erste Algorithmus, der in den Sinn kommt, ist der Goertzel-Algorithmus . Dieser Algorithmus geht normalerweise davon aus, dass die interessierende Frequenz ein ganzzahliges Vielfaches der Grundfrequenz ist. In diesem Artikel wird jedoch der (verallgemeinerte) Algorithmus auf den Fall angewendet, an dem Sie interessiert sind.


Ein weiteres Problem ist, dass das Signalmodell falsch ist. Es nutzt 2*%pi*(1:siglen)*(Fc/siglen). Es sollte verwendet werden, 2*%pi*(0:siglen-1)*(Fc/siglen)damit die Phase korrekt ausgegeben wird.

Ich denke auch, dass es ein Problem mit der Frequenz gibt Fc=21.3, die sehr niedrig ist. Niedrigfrequente reelle Signale neigen dazu, eine Vorspannung aufzuweisen, wenn es um Probleme bei der Phasen- / Frequenzschätzung geht.

Ich habe auch eine Grobgittersuche nach der Phasenschätzung durchgeführt, die dieselbe Antwort liefert wie der Goertzel-Algorithmus.

Unten sehen Sie eine grafische Darstellung, die die Abweichung in beiden Schätzungen (Goertzel: blau, Coarse: rot) für zwei verschiedene Frequenzen zeigt: Fc=21.3(durchgehend) und Fc=210.3(gestrichelt). Wie Sie sehen können, ist die Vorspannung für die höhere Frequenz viel geringer.

Die Kurve Achse ist die Anfangsphase, die sich von 0 auf 2 π ändert .x2π

Bildbeschreibung hier eingeben

Peter K.
quelle
Ich habe gerade den Code für den Goerzel-Algorithmus getestet, der auf dem Papier basiert. Mit dem ausgegebenen DTFT-Wert kann die Spitze sehr genau erhalten werden. Es gibt jedoch einen Skalierungsfaktor von genau 1000. Wenn also der ursprüngliche Höchstwert 1.234 beträgt, ist er nach Goerzel 1234. Weiß jemand, woher dies kommen könnte?
IR8N6I
Hat in der Zwischenzeit nachgeforscht. Wahrscheinlich hat es mit der Amplitudenskalierung zu tun: Skalierung der Zeitdomänenamplitude = Frequenzdomänenkoeffizient * 2 / N, wobei N die Länge des Signals ist. Ist diese Annahme richtig?
IR8N6I
Hallo! Ich habe gerade herausgefunden, dass mit dem Goertzel-Algorithmus die Amplitude beim resultierenden komplexen Koeffizienten sehr genau ist, aber die Phase ist völlig falsch. Hat jemand eine Idee, woher das kommen könnte? Mit "Phase" meine ich die Phasenverzögerung, die in der Grundwelle des ursprünglichen Signals angegeben ist.
IR8N6I
1
sin(ω0t+ϕ)j2[ejϕδ~(ω+ω0+2πk)e+jϕδ~(ωω0+2πk)]π/2
4

Wenn Sie mehrere benachbarte FFT-Bins verwenden möchten, nicht nur 2, kann die fensterbasierte Sinc-Interpolation zwischen den Ergebnissen der komplexen Bins abhängig von der Breite des Fensters eine sehr genaue Schätzung ergeben.

Windowed-Sinc-Interpolation ist häufig in hochwertigen Audio-Upsamplern zu finden, sodass Artikel zu diesem Thema über geeignete Interpolationsformeln mit Fehleranalyse verfügen.

hotpaw2
quelle
Danke für den Kommentar. Ich werde diesen Ansatz auch ausprobieren.
IR8N6I
4

sin(πx)(πx)

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

[2] K. Dressler, „Sinusförmige Extraktion unter Verwendung einer effizienten Implementierung einer FFT mit mehreren Auflösungen“, in Proc. 9. Int. Conf. über digitale Audioeffekte (DAFx-06), Montreal, Kanada, September 2006, S. 247–252.

ederwander
quelle
Hallo! Vielen Dank für alle Ihre Kommentare. Ich habe meinen Code (siehe unten) erweitert, um Goertzel-Filter mit parabolischer Spitzeninterpolation zu kombinieren, um die Phase zu erhalten. Die Ergebnisse sind jedoch immer noch nicht genau (+ - 3-4 Grad). Ist das so nah wie möglich oder gibt es Fehler beim Verstehen oder Codieren?
IR8N6I
3

Vor ein paar Jahren hatte ich große Schwierigkeiten mit diesem Problem.

Ich habe diese Frage gestellt:

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

Am Ende habe ich die Berechnungen von Grund auf neu durchgeführt und eine Antwort auf meine eigene Frage gepostet.

Ich bin überrascht, dass ich im Internet keine ähnliche Ausstellung finden konnte.

Ich werde die Antwort hier noch einmal posten. Beachten Sie, dass der Code für ein Szenario entwickelt wurde, in dem ich mein FFT-Fenster um das Vierfache überlappe.

π


Dieses Puzzle benötigt zwei Schlüssel, um es freizuschalten.

Grafik 3.3:

Bildbeschreibung hier eingeben

Grafik 3.4:

Bildbeschreibung hier eingeben

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
quelle
Sie interpolieren die Frequenz, während das OP die Frequenz kennt und die Amplitude interpolieren möchte.
15.
2

Mit diesem Python-Code erhalten Sie ein sehr genaues Ergebnis (ich habe es für viele Noten verwendet und Fehler von weniger als 0,01% des Halbtons erhalten) mit parabolischer Interpolation (Methode erfolgreich von McAulay Quatieri, Serra usw. in Harmonic + Residual verwendet) Trenntechniken)

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

Die Frequenzen, mit denen Sie es zu tun haben (21,3 Hz, abgetastet bei 8 kHz), sind sehr niedrig. Da es sich um reelle Signale handelt, weisen sie für ** jede ** Frequenz eine Verzerrung in der Phasenschätzung auf.

Dieses Bild zeigt eine grafische Darstellung der Vorspannung ( phase_est - phase_orig) für Fc = 210.3;(in Rot) gegenüber der Vorspannung für Fc = 21.3;. Wie Sie sehen, ist der Versatz für den 21.3Fall viel wichtiger .

Eine weitere Option ist die Reduzierung Ihrer Abtastrate. Die grüne Kurve zeigt die Abweichung Fs = 800von 8000.

Bildbeschreibung hier eingeben

IR8n6i
quelle
1
Danke für das Update! Siehe meine Handlung; Ich denke immer noch, dass jeder Phasenschätzer einen Bias für diese niedrige Frequenz haben wird. Eine Möglichkeit, dies zu umgehen, besteht darin, die bekannte Frequenz (falls bekannt!) Zu verwenden, um die Phasenschätzungsvorspannung durch eine Nachschlagetabelle zu korrigieren. Aber Sie müssen vorsichtig sein: Die Vorspannung ändert sich mit der Frequenz. Eine andere Möglichkeit besteht darin, die Abtastrate zu verringern.
Peter K.
1
Danke dir auch! Wenn Sie jedoch Fs = 8000 Hz und Fc = 210 anstelle von 210,3 verwenden, sieht die Vorspannung noch schlechter aus. Irgendeine Idee, woher das kommen könnte?
IR8N6I
1
Erk! Keine Ahnung. FWIW, wird der Geortzel Schätzer Probleme nicht haben: goertzel = atan(imag(y(2)),real(y(2)))*180/%pi + 90;. :-) Werde noch ein bisschen graben. Beobachten Sie diesen Raum.
Peter K.
1
Die parabolische Interpolation entspricht nicht Ihren Vorstellungen. Insbesondere wenn Sie Ihre Berechnung von pdurch ersetzen, erhalten 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;Sie VIEL bessere Antworten - auch für Fc=210. Ich bin mir nicht sicher, ob die aktuelle Version von pIhnen etwas Sinnvolles bringt. Die Interpolationsformel dient zur Interpolation der AMPLITUDE einer Parabel, pinterpoliert aber die gerade ... ungerade Phase .
Peter K.
1
Das ist alles in Ordnung, außer dass die Peak-Position ( p = (yp1 - ym1)/(2*(2*y0 - yp1 - ym1))) manchmal falsch ist, wenn Sie PHASES anstelle von Amplituden verwenden. Dies liegt daran, dass die Phasen um die Grenze von +/- 180 Grad springen können . Alles, was benötigt wird, um es für die Phase zu reparieren, ist, diese Linie zu meiner p2obigen Berechnung zu ändern .
Peter K.