Czy istnieje algorytm obliczający fazę dla pojedynczej częstotliwości?

Jeśli masz funkcję i odwołujesz się do sin wave co by to był szybki algorytm do obliczenia ?sin ( ω x ) $f(t)=A \cdot \sin(\omega t+\phi)$$\sin(\omega x)$$\phi$

Szukałem na Goertzela algorytmu, ale nie wydaje się do czynienia z fazą?

Użyj DFT na określonej częstotliwości. Następnie oblicz amplitudę i fazę z części rzeczywistych / imagowych. Daje ci fazę odniesioną do początku czasu próbkowania.

W „normalnym” FFT (lub DFT obliczonym dla wszystkich N harmonicznych) zwykle obliczasz częstotliwość za pomocą f = k * (sample_rate) / N, gdzie k jest liczbą całkowitą. Chociaż może to wydawać się świętokradztwo (szczególnie dla członków Kościoła z Całkowitą Liczbą Całkowitą), możesz faktycznie używać wartości k nie będących liczbami całkowitymi podczas wykonywania pojedynczego DFT.

Załóżmy na przykład, że wygenerowałeś (lub uzyskałeś) N = 256 punktów fali sinusoidalnej o częstotliwości 27 Hz. (powiedzmy, sample_rate = 200). Twoje „normalne” częstotliwości dla 256-punktowego FFT (lub punktu N DFT) odpowiadają: f = k * (próbka_rate) / N = k * (200) / 256, gdzie k jest liczbą całkowitą. Ale liczba całkowita „k” nie będąca liczbą całkowitą wynosząca 34,56 odpowiadałaby częstotliwości 27 Hz, przy użyciu parametrów wymienionych powyżej. To tak, jakby stworzyć „bin” DFT, który jest dokładnie wyśrodkowany na częstotliwości zainteresowania (27 Hz). Niektóre kody C ++ (kompilator DevC ++) mogą wyglądać następująco:

#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: Mam nadzieję, że powyższe powyższe dobrze przekłada się na przepełnienie stosu - niektóre z nich mogą się owijać)

Wynikiem powyższego jest faza -twopi / 4, jak pokazano w wygenerowanych punktach rzeczywistych (a wzmacniacz jest podwojony, aby odzwierciedlić częstotliwość dodatnią / ujemną).

Należy zwrócić uwagę na kilka rzeczy - używam cosinusa do generowania fali testowej i interpretowania wyników - musisz być ostrożny - faza odnosi się do czasu = 0, czyli od momentu rozpoczęcia próbkowania (tj. Kiedy zebrałeś r [0] ), a cosinus jest poprawną interpretacją).

Powyższy kod nie jest elegancki ani wydajny (np .: użyj tabel przeglądowych dla wartości sin / cos itp.).

Twoje wyniki będą dokładniejsze, gdy użyjesz większego N, i jest trochę błędu ze względu na fakt, że częstotliwość próbkowania i N powyżej nie są wielokrotnościami.

Oczywiście, jeśli chcesz zmienić częstotliwość próbkowania, N lub f, musisz zmienić kod i wartość k. Możesz umieścić bin DFT w dowolnym miejscu na linii ciągłej częstotliwości - po prostu upewnij się, że używasz wartości k, która odpowiada częstotliwości zainteresowania.

Problem można sformułować jako (nieliniowy) problem najmniejszych kwadratów:

gdzie to funkcja celu, którą należy zminimalizować w odniesieniu do .$F(\phi)$$\phi$

Pochodna jest bardzo prosta:

Powyższa funkcja cel może zostać zminimalizowany iteracyjnie stosując metody największego spadku (aproksymacja pierwszego rzędu), metoda Newtona , metody Gaussa-Newtona lub metody Levenberga Marquardt (drugi celu aproksymacji - muszą być dostarczane w nich).$F''(\phi)$

Oczywiście powyższa funkcja celu ma wiele minimów z powodu okresowości, dlatego też można dodać pewien okres karny w celu rozróżnienia innych minimów (na przykład dodając do równania modelu). Ale myślę, że optymalizacja będzie właśnie zbiegają się z dokładnością do minimów i można zaktualizować wynik odejmowania . 2 π k , k ∈ $\phi^{2}$$2\pi k, k\in N$

Istnieje kilka różnych sformułowań algorytmu Goertzela. Te, które zapewniają 2 zmienne stanu (ortogonalne lub bliskie) lub złożone zmienne stanu, jako możliwe dane wyjściowe, często można wykorzystać do obliczenia lub oszacowania fazy w odniesieniu do jakiegoś punktu w oknie Goertzela, takiego jak środek. Te, które zapewniają tylko jeden wynik skalarny, zwykle nie mogą.

Musisz także wiedzieć, gdzie znajduje się twoje okno Goertzela w stosunku do osi czasu.

Jeśli twój sygnał nie jest dokładnie liczbą całkowitą okresową w oknie Goertzela, oszacowanie fazy wokół punktu odniesienia na środku okna może być dokładniejsze niż odniesienie fazy do początku lub końca.

Pełna FFT to przesada, jeśli znasz częstotliwość swojego sygnału. Ponadto Goertzel można dostroić do częstotliwości nieokresowej na długości FFT, podczas gdy FFT będzie wymagało dodatkowej interpolacji lub zerowania dla częstotliwości nieokresowych w oknie.

Złożony Goertzel jest równoważny 1 binowi DFT, który wykorzystuje wznowę dla wektorów cosinus i sinus lub czynników skręcających FFT.

Jeśli twoje sygnały są wolne od szumów, możesz zidentyfikować przejścia przez zero w obu i określić częstotliwość i fazę względną.

Zależy to od tego, jaka jest twoja definicja „szybkiego”, jak dokładna jest twoja ocena, czy chcesz lub fazę względem twoich próbek, i ile hałasu jest na twojej funkcji i fali sinusoidalnej.$\phi$

Jednym ze sposobów jest zrobienie FFT z i po prostu spójrz na bin najbliższy . $f(t)$$\omega$ Będzie to jednak zależeć od tego, czy jest zbliżone do środkowej częstotliwości bin.$\omega$

Więc:

• Co rozumiesz przez „szybki”?
• Jak dokładny jest szacunek?
• Czy chcesz (faza względem odniesienia) czy faza względem początku próbkowania? Czy to ma znaczenie?$\phi$
• Jaki jest poziom hałasu dla każdego sygnału?

PS: Zakładam, że miałeś na myśli , a nie .$f(t)=A\sin(ωt+ϕ)$$f({\Huge t})=A\sin(ω{\Huge x}+ϕ)$

Punkt początkowy:
1) pomnóż falę sinusoidalną sygnału i odniesienia: = A⋅sin (ωt + ϕ) ⋅sin (ωt) = 0,5⋅A⋅ (cos (ϕ) - cos (2⋅ωt + ϕ) ) 2) znajdź całkę z okresu : 3) możesz obliczyć :
$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)$

Pomyśl:
jak zmierzyć A?
jak określić w przedziale ? (pomyśl o „referencyjnej fali cos ”)$\phi$$0..(2 \cdot \pi)$

W przypadku sygnału dyskretnego zmień całkę na sumę i ostrożnie wybierz T!

Możesz także to zrobić (w notacji numerycznej):

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

gdzie sygnał jest sygnałem przesuniętym w fazie, cos i sin są sygnałami odniesienia, a ty generujesz aproksymację całki w określonym czasie poprzez zsumowanie dwóch produktów.

Jest to poprawa sugestii @Kevin McGee, aby zastosować DFT o pojedynczej częstotliwości z ułamkowym indeksem bin. Algorytm Kevina nie daje świetnych rezultatów: podczas gdy przy półkach i całych przedziałach jest bardzo dokładny, również blisko dziur i połówek, jest również całkiem dobry, ale w przeciwnym razie błąd może wynosić 5%, co prawdopodobnie nie jest do zaakceptowania dla większości zadań .

Proponuję ulepszyć algorytm Kevina, dostosowując , tj. Długość okna DFT, aby zbliżyło się do całości, jak to możliwe. Działa to, ponieważ w przeciwieństwie do FFT, DFT nie wymaga, aby było potęgą 2.$N$$k$$N$

Poniższy kod jest w Swift, ale powinien być intuicyjnie zrozumiały:

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