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

17

Jeśli masz funkcję i odwołujesz się do sin wave co by to był szybki algorytm do obliczenia ?sin ( ω x ) ϕf(t)=Asin(ωt+ϕ)sin(ωx)ϕ

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

SamFisher83
źródło

Odpowiedzi:

5

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.

Kevin McGee
źródło
To podejście można ulepszyć, dostosowując N, aby zbliżyć k do całości. Opublikowałem osobną odpowiedź, która poprawia dokładność tego algorytmu.
mojuba,
10

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

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

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

Pochodna jest bardzo prosta:

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

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

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

Libor
źródło
Nie sądzę, że musisz karać z powodu okresowości nie? Możesz wziąć dowolne minima w przestrzeni fazowej, z którą się zbiega i zrobić modulu , nie? 2π
Spacey,
@Mohammad Tak, ale niektóre techniki optymalizacji mogą wykorzystywać wiele punktów początkowych, które powinny zbiegać się do tej samej wartości lub przyjmować funkcję wypukłą z pojedynczym globalnym minimalizatorem, który można dobrze przybliżać za pomocą kwadratu. Inną korzyścią jest to, że kończymy z takim samym wynikiem dla dowolnego punktu początkowego . ϕ0
Libor,
Ciekawy. Czy mogę zaprosić Cię również do odpowiedzi na to powiązane pytanie ? :-)
Spacey
@Mohammad OK, przyczyniłem się trochę tam :)
Libor,
Dokąd idzie funkcja fi (w)? fi (w) nie jest stałą, więc kiedy weźmiesz pochodną niestałej, jak staje się zero?
SamFisher83,
5

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.

hotpaw2
źródło
Czy oszacowanie fazy nigdzie w oknie nie jest dokładnie tej samej dokładności, ponieważ wystarczy dodać do oszacowania fazy na początku okna, aby obliczyć oszacowanie fazy dla próbki w oknie ( będący początkiem okna)? k k = 0ωkkk=0
Olli Niemitalo
Nie, ponieważ dodanie wk powoduje inną fazę na końcu okna niż na początku dla sinusoidy nieokreślonej okresowo w aperturze. Ale 1-bin DFT oblicza pojedynczą fazę kołową w tym samym punkcie. Zatem wszystkie 3 wartości będą różne. Ale faza środkowa jest zawsze związana ze stosunkiem funkcji nieparzystej / parzystej, bez względu na f0.
hotpaw2
Próbuję, ale nie rozumiem.
Olli Niemitalo
Użyj cosinusa (faza zero przy k = 0), nieznacznie podkręć częstotliwość (o małą liczbę irracjonalną, ale bez zmiany fazy przy k = 0). DFT informuje, że faza uległa zmianie! Spróbuj tego samego z cosinusem dokładnie wyśrodkowanym na k = N / 2. Bez zmian przy k = N / 2 dla dowolnego df. To samo dotyczy grzechu lub dowolnej mieszanki. Centrowanie punktu odniesienia fazy pokazuje mniej zmian w zmierzonej fazie wraz ze zmianami w f0. np. błąd częstotliwości nie przyczynia się do zwiększenia błędów pomiaru fazy.
hotpaw2
1
Tak, błąd oszacowania fazy mniejszy w środku okna ma sens, jeśli sinusoida i filtr Goertzela mają różne częstotliwości. W takim przypadku oszacowanie fazy, powiedziane na końcu okna, jest tendencyjne przez stałą, która jest iloczynem odległości między środkiem i końcem okna oraz różnicy między częstotliwościami filtra sinusoidalnego i Goertzela. Odejmowanie tego odchylenia daje ten sam błąd wielkości, co dla oszacowania środkowego, ale wymaga znajomości częstotliwości sinusoidy.
Olli Niemitalo
4

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

Juancho
źródło
3

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

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

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?ϕ
  • Jaki jest poziom hałasu dla każdego sygnału?

PS: Zakładam, że miałeś na myśli , a nie .f(t)=Asin(ωt+ϕ)f(t)=Asin(ωx+ϕ)

Peter K.
źródło
2

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=π/ω
I(ϕ)=0TF(t)dt =0.5Acos(ϕ)T
ϕ
cos(ϕ)=I(t)/(0.5AT)

Pomyśl:
jak zmierzyć A?
jak określić w przedziale ? (pomyśl o „referencyjnej fali cos ”)ϕ0..(2π)

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

SergV
źródło
1

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.

M529
źródło
0

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

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))
mojuba
źródło
FFT to po prostu sposób na skuteczne obliczenie DFT. W nowoczesnych bibliotekach nie ma już mocy dwóch ograniczeń. Jeśli potrzebujesz tylko jednej lub dwóch wartości bin, lepiej obliczyć je bezpośrednio, tak jak zrobiłeś. Dla pojedynczego czystego tonu (rzeczywistego lub złożonego) potrzebne są tylko dwie wartości bin, aby dokładnie obliczyć częstotliwość, fazę i amplitudę. Zobacz dsprelated.com/showarticle/1284.php . Matematyka jest dość skomplikowana, ale istnieją linki do artykułów, w których wyjaśniono pochodne. Algebra liniowa jest warunkiem prawdziwego zrozumienia.
Cedron Dawg,