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

17

Wenn Sie eine Funktion f(t)=Asin(ωt+ϕ) und eine Referenz-Sin-Welle sin(ωx) was wäre ein schneller Algorithmus zur Berechnung von ϕ ?

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

SamFisher83
quelle

Antworten:

5

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.

Kevin McGee
quelle
Dieser Ansatz kann verbessert werden, indem N angepasst wird, um k näher an ein Ganzes heranzuführen. Ich habe eine separate Antwort veröffentlicht, die die Genauigkeit dieses Algorithmus verbessert.
Mojuba
10

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

F(ϕ)=12i=1n[Asin(ωi+ϕ)fi(ω)]2

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

Die Ableitung ist sehr einfach:

F(ϕ)=i=1nAcos(ωi+ϕ)[Asin(ωi+ϕ)fi(ω)]

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(ϕ)

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 .ϕ22πk,kN

Libor
quelle
Ich glaube nicht, dass Sie wegen der Periodizität bestraft werden müssen, nein? Sie können einfach alle Minima im Phasenraum nehmen, zu denen es konvergiert, und ein Modul ausführen , nein? 2π
Spacey
@Mohammad Ja, aber einige Optimierungstechniken verwenden möglicherweise mehrere Startpunkte, die auf denselben Wert konvergieren sollten, oder nehmen eine konvexe Funktion mit einem einzelnen globalen Minimierer an, der mit einem Quadrat gut angenähert werden kann. Der andere Vorteil ist, dass wir für jeden Startpunkt mit demselben Ergebnis enden . ϕ0
Libor
Interessant. Darf ich Sie einladen, sich auch mit dieser verwandten Frage zu befassen? :-)
Spacey
@Mohammad OK, ich habe dort ein wenig beigetragen :)
Libor
Wohin geht die Funktion fi (w)? fi (w) ist keine Konstante. Wenn Sie also eine Ableitung einer Nichtkonstante nehmen, wie wird sie dann Null?
SamFisher83
5

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.

hotpaw2
quelle
ωkkk=0
Nein, da das Hinzufügen von wk am Ende des Fensters zu einer anderen Phase führt als zu Beginn für eine Sinuskurve mit nicht ganzzahliger periodischer Apertur. Eine 1-Bin-DFT berechnet jedoch an demselben Punkt eine einzelne Kreisphase. Somit sind die 3 Werte alle unterschiedlich. Die Mittenphase hängt jedoch immer vom Verhältnis der ungeraden / geraden Funktion ab, unabhängig von f0.
hotpaw2
Ich versuche es, aber das verstehe ich nicht.
Olli Niemitalo
Verwenden Sie einen Kosinus (Phase Null bei k = 0), passen Sie die Frequenz leicht an (um eine winzige irrationale Zahl, ohne jedoch die Phase bei k = 0 zu ändern). Eine DFT meldet, dass sich die Phase geändert hat! Versuchen Sie dasselbe mit einem Kosinus, der genau bei k = N / 2 zentriert ist. Keine Änderung bei k = N / 2 für irgendeinen df. Gleiches gilt für Sünde oder irgendeine Mischung. Das Zentrieren des Phasenreferenzpunkts zeigt weniger Änderungen in der gemessenen Phase mit Änderungen in f0. zB trägt ein Frequenzfehler nicht zu erhöhten Phasenmessfehlern bei.
hotpaw2
1
Ja, der Phasenschätzungsfehler in der Mitte des Fensters ist sinnvoll, wenn Sinus und Goertzel-Filter unterschiedliche Frequenzen haben. In diesem Fall wird die Phasenschätzung beispielsweise am Ende des Fensters durch eine Konstante vorgespannt, die das Produkt aus dem Abstand zwischen der Mitte und dem Ende des Fensters und der Differenz zwischen den Sinus- und Goertzel-Filterfrequenzen ist. Das Subtrahieren dieser Vorspannung ergibt den gleichen Größenfehler wie für die Mittenschätzung, erfordert jedoch die Kenntnis der Frequenz der Sinuskurve.
Olli Niemitalo
4

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

Juancho
quelle
3

ϕ

f(t)ωω

So:

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

f(t)=Asin(ωt+ϕ)f(t)=Asin(ωx+ϕ)

Peter K.
quelle
2



F(t)
T=π/ω
I(ϕ)=0TF(t)dt =0.5Acos(ϕ)T
ϕ
cos(ϕ)=I(t)/(0.5AT)



ϕ0..(2π)

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

SergV
quelle
1

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.

M529
quelle
0

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 .

NkN

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))
Mojuba
quelle
Die FFT ist einfach eine Möglichkeit, eine DFT effizient zu berechnen. Mit modernen Bibliotheken ist die Kraft von zwei Einschränkungen nicht mehr vorhanden. Wenn Sie nur ein oder zwei Bin-Werte benötigen, ist es besser, diese direkt wie Sie zu berechnen. Für einen einzelnen reinen Ton (real oder komplex) werden nur zwei Bin-Werte benötigt, um Frequenz, Phase und Amplitude genau zu berechnen. Siehe dsprelated.com/showarticle/1284.php . Die Mathematik ist ziemlich ausgefeilt, aber es gibt Links zu den Artikeln, in denen die Ableitungen erklärt werden. Die lineare Algebra ist eine Voraussetzung für ein echtes Verständnis.
Cedron Dawg