Gibt es einen Algorithmus zur Berechnung der Phase für eine einzelne Frequenz?

Wenn Sie eine Funktion $f(t)=A \cdot \sin(\omega t+\phi)$ und eine Referenz-Sin-Welle $\sin(\omega x)$ was wäre ein schneller Algorithmus zur Berechnung von $\phi$ ?

Ich habe mir den Goertzel- Algorithmus angesehen, aber er scheint sich nicht mit der Phase zu befassen.

Verwenden Sie eine DFT mit der spezifischen Frequenz. Berechnen Sie dann Amplitude und Phase aus den Real / Imag-Teilen. Sie erhalten die Phase, die auf den Beginn der Abtastzeit bezogen ist.

In einer 'normalen' FFT (oder einer für alle N Harmonischen berechneten DFT) berechnen Sie normalerweise die Frequenz mit f = k * (sample_rate) / N, wobei k eine ganze Zahl ist. Obwohl es sakrilegisch erscheint (insbesondere für Mitglieder der Church of the Wholly Integer), können Sie bei einer einzelnen DFT tatsächlich nicht ganzzahlige Werte von k verwenden.

Angenommen, Sie haben N = 256 Punkte einer Sinuswelle von 27 Hz erzeugt (oder erhalten). (Nehmen wir an, sample_rate = 200). Ihre 'normalen' Frequenzen für eine 256-Punkt-FFT (oder N-Punkt-DFT) würden entsprechen: f = k * (Abtastrate) / N = k * (200) / 256, wobei k eine ganze Zahl ist. Ein nicht ganzzahliges 'k' von 34,56 würde jedoch einer Frequenz von 27 Hz unter Verwendung der oben aufgeführten Parameter entsprechen. Es ist so, als würde man einen DFT-Bin erstellen, der genau auf der interessierenden Frequenz (27 Hz) zentriert ist. Einige C ++ - Codes (DevC ++ - Compiler) sehen möglicherweise wie folgt aus:

#include <cstdio>
#include <cstdlib>
#include <iostream>
#include <cmath>
using namespace std;

// arguments in main needed for Dev-C++ I/O
int main (int nNumberofArgs, char* pszArgs[ ] ) {
const long N = 256 ;
double sample_rate = 200., amp, phase, t, C, S, twopi = 6.2831853071795865; 
double  r[N] = {0.}, i[N] = {0.}, R = 0., I = 0. ;
long n ;

// k need not be integer
double k = 34.56;

// generate real points
for (n = 0; n < N; n++) {
    t =  n/sample_rate;
    r[n] = 10.*cos(twopi*27.*t - twopi/4.);
}  // end for

// compute one DFT
for (n = 0; n < N; n++) {
    C = cos(twopi*n*k/N); S = sin(twopi*n*k/N);
    R = R + r[n]*C + i[n]*S;
    I = I + i[n]*C - r[n]*S;
} // end for

cout<<"\n\ndft results for N = " << N << "\n";
cout<<"\nindex k     real          imaginary       amplitude         phase\n";

amp = 2*sqrt( (R/N)*(R/N) + (I/N)*(I/N) ) ;
phase = atan2( I, R ) ;
// printed R and I are scaled
printf("%4.2f\t%11.8f\t%11.8f\t%11.8f\t%11.8f\n",k,R/N,I/N,amp,phase);

cout << "\n\n";
system ("PAUSE");
return 0;
} // end main

//**** end program

(PS: Ich hoffe, dass das oben Gesagte gut zu einem Stapelüberlauf führt - ein Teil davon könnte sich umwickeln)

Das Ergebnis des Obigen ist eine Phase von -twopi / 4, wie in den erzeugten realen Punkten gezeigt (und der Verstärker wird verdoppelt, um die pos / neg-Frequenz widerzuspiegeln).

Ein paar Dinge zu beachten - ich benutze Cosinus, um die Testwellenform zu generieren und die Ergebnisse zu interpretieren - Sie müssen vorsichtig sein - Phase bezieht sich auf Zeit = 0, wenn Sie mit dem Sampling begonnen haben (dh wenn Sie r [0] gesammelt haben) ) und Cosinus ist die richtige Interpretation).

Der obige Code ist weder elegant noch effizient (z. B. Verwenden Sie Nachschlagetabellen für die sin / cos-Werte usw.).

Ihre Ergebnisse werden genauer, wenn Sie ein größeres N verwenden, und es liegt ein kleiner Fehler vor, da die Abtastrate und das obige N kein Vielfaches voneinander sind.

Wenn Sie Ihre Abtastrate N oder f ändern möchten, müssen Sie natürlich den Code und den Wert von k ändern. Sie können einen DFT-Bin an einer beliebigen Stelle auf der durchgehenden Frequenzlinie ablegen. Stellen Sie lediglich sicher, dass Sie einen Wert von k verwenden, der der interessierenden Frequenz entspricht.

Das Problem kann als (nichtlineares) Problem der kleinsten Quadrate formuliert werden:

wobei die Zielfunktion ist, die in Bezug auf ϕ zu minimieren ist .$F(\phi)$$\phi$

Die Ableitung ist sehr einfach:

Die obige Zielfunktion kann iterativ unter Verwendung der Gradientenabstiegsmethode (Näherung erster Ordnung), der Newton-Methode , der Gauß-Newton-Methode oder der Levenberg-Marquardt-Methode (Näherung zweiter Ordnung - muss in diesen bereitgestellt werden) minimiert werden .$F''(\phi)$

Offensichtlich hat die obige Zielfunktion aufgrund der Periodizität mehrere Minima, daher kann ein Strafbegriff hinzugefügt werden, um andere Minima zu unterscheiden (zum Beispiel Hinzufügen von zur Modellgleichung). Ich denke jedoch, dass die Optimierung nur auf die nächsten Minima konvergiert und Sie das Ergebnis aktualisieren können, indem Sie 2 π k , k ∈ N abziehen .$\phi^{2}$$2\pi k, k\in N$

Es gibt verschiedene Formulierungen des Goertzel-Algorithmus. Diejenigen, die 2 Zustandsvariablen (orthogonal oder nahe beieinander) oder eine komplexe Zustandsvariable als mögliche Ausgaben bereitstellen, können häufig verwendet werden, um die Phase unter Bezugnahme auf einen Punkt im Goertzel-Fenster, wie z. B. die Mitte, zu berechnen oder zu schätzen. Diejenigen, die nur eine einzige skalare Ausgabe liefern, können dies normalerweise nicht.

Sie müssen auch wissen, wo sich Ihr Goertzel-Fenster in Bezug auf Ihre Zeitachse befindet.

Wenn Ihr Signal in Ihrem Goertzel-Fenster nicht genau ganzzahlig ist, ist die Phasenschätzung um einen Referenzpunkt in der Mitte des Fensters möglicherweise genauer als die Referenzierungsphase auf den Anfang oder das Ende.

Eine vollständige FFT ist übertrieben, wenn Sie die Frequenz Ihres Signals kennen. Außerdem kann ein Goertzel auf eine Frequenz abgestimmt werden, die in der FFT-Länge nicht periodisch ist, während eine FFT für nicht periodische Frequenzen im Fenster zusätzliche Interpolation oder Null-Auffüllung benötigt.

Ein komplexes Goertzel entspricht 1 Bin einer DFT, die eine Wiederholung für die Cosinus- und Sinusbasisvektoren oder FFT-Twiddle-Faktoren verwendet.

Wenn Ihre Signale rauschfrei sind, können Sie Nulldurchgänge in beiden identifizieren und die Frequenz und die relative Phase bestimmen.

$\phi$

$f(t)$$\omega$$\omega$

So:

• Was meinst du mit "schnell"?
• Wie genau benötigen Sie die Schätzung?
• $\phi$
• Wie hoch ist der Rauschpegel bei jedem Signal?

$f(t)=A\sin(ωt+ϕ)$$f({\Huge t})=A\sin(ω{\Huge x}+ϕ)$

$F(t)$
$T= \pi /\omega$
$I(\phi) = \int_0^T F(t)dt\ = 0.5⋅A⋅cos(ϕ) \cdot T$
$\phi$
$cos(\phi) = I(t)/(0.5 \cdot A \cdot T)$

$\phi$$0..(2 \cdot \pi)$

Für diskrete Signale ändern Sie das Integral in Summe und wählen Sie sorgfältig T!

Sie können dies auch tun (in Numpy-Notation):

np.arctan( (signal*cos).sum() / (signal*sin).sum() ))

Wenn das Signal Ihr phasenverschobenes Signal ist, sind cos und sin die Referenzsignale, und Sie erzeugen eine Annäherung eines Integrals über eine bestimmte Zeit durch Summieren über die beiden Produkte.

Dies ist eine Verbesserung gegenüber @ Kevin McGees Vorschlag, eine Einzelfrequenz-DFT mit einem gebrochenen Bin-Index zu verwenden. Kevins Algorithmus liefert keine großartigen Ergebnisse: Während er bei halben und ganzen Behältern sehr präzise ist, ist er auch in der Nähe der ganzen und halben Hälften ziemlich gut, aber ansonsten kann der Fehler innerhalb von 5% liegen, was für die meisten Aufgaben wahrscheinlich nicht akzeptabel ist .

$N$$k$$N$

Der folgende Code ist in Swift, sollte aber intuitiv klar sein:

let f = 27.0 // frequency of the sinusoid we are going to generate
let S = 200.0 // sampling rate
let Nmax = 512 // max DFT window length
let twopi = 2 * Double.pi

// First, calculate k for Nmax, and then round it
var k = round(f * Double(Nmax) / S)

// The magic part: recalculate N to make k as close to whole as possible
// We also need to recalculate k once again due to rounding of N. This is important.
let N = Int(k * S / f)
k = f * Double(N) / S

// Generate the sinusoid
var r: [Double] = []
for i in 0..<N {
    let t = Double(i) / S
    r.append(sin(twopi * f * t))
}

// Compute single-frequency DFT
var R = 0.0, I = 0.0
let twopikn = twopi * k / Double(N)
for i in 0..<N {
    let x = Double(i) * twopikn
    R += r[i] * cos(x)
    I += r[i] * sin(x)
}
R /= Double(N)
I /= Double(N)

let amp = 2 * sqrt(R * R + I * I)
let phase = atan2(I, R) / twopi

print(String(format: "k = %.2f    R = %.8f    I = %.8f    A = %.8f    φ/2π = %.8f", k, R, I, amp, phase))