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

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.

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 .$x$$2\pi$

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.

[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.

Eine Methode besteht darin, das Maximum zu finden und eine Parabel darüber anzupassen und dann das Maximum des Parabolas als Frequenz- und Größenschätzung zu verwenden. Hier können Sie alles darüber lesen: https://ccrma.stanford.edu/~jos/sasp/Sinusoidal_Peak_Interpolation.html

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.

• Der erste Schlüssel ist zu verstehen, wie das Überlappen des FFT-Fensters eine Rotation in der Bin-Phase einführt.

• Der zweite Schlüssel stammt aus Grafik 3.3 & 3.4 hier (Dank an Stephan Bernsee für die Erlaubnis, die Bilder hier zu kopieren).

Grafik 3.3:

Grafik 3.4:

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;
}

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'

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.