Apakah ada algoritma untuk menghitung fase untuk frekuensi tunggal?

Jika Anda memiliki fungsi , dan referensi dosa gelombang dosa ( ω x ) apa yang akan menjadi algoritma cepat untuk menghitung ϕ ?$f(t)=A \cdot \sin(\omega t+\phi)$$\sin(\omega x)$$\phi$

Saya melihat algoritma Goertzel , tetapi sepertinya tidak berurusan dengan fase?

Gunakan DFT pada frekuensi tertentu. Kemudian hitung amplitudo dan fase dari bagian real / imag. Ini memberi Anda fase yang dirujuk ke awal waktu pengambilan sampel.

Dalam FFT 'normal' (atau DFT dihitung untuk semua harmonik N), Anda biasanya menghitung frekuensi dengan f = k * (sample_rate) / N, di mana k adalah bilangan bulat. Meskipun mungkin tampak asusila (terutama untuk anggota Gereja Wholly Integer), Anda sebenarnya dapat menggunakan nilai-nilai non-integer k ketika melakukan DFT tunggal.

Sebagai contoh, misalkan Anda telah menghasilkan (atau memperoleh) N = 256 poin dari gelombang sinus 27 Hz. (katakanlah, sample_rate = 200). Frekuensi 'normal' Anda untuk 256 titik FFT (atau N titik DFT) akan sesuai dengan: f = k * (sample_rate) / N = k * (200) / 256, di mana k adalah bilangan bulat. Tetapi 'k' non-integer dari 34,56 akan sesuai dengan frekuensi 27 Hz., Menggunakan parameter yang tercantum di atas. Ini seperti membuat 'bin' DFT yang persisnya berpusat pada frekuensi yang diinginkan (27 Hz.). Beberapa kode C ++ (kompilator DevC ++) mungkin terlihat sebagai berikut:

#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: Saya harap di atas diterjemahkan dengan baik untuk stackoverflow - beberapa mungkin membungkus)

Hasil di atas adalah fase -twopi / 4, seperti yang ditunjukkan pada titik nyata yang dihasilkan (dan amp digandakan untuk mencerminkan frekuensi pos / neg).

Beberapa hal yang perlu diperhatikan - saya menggunakan cosinus untuk menghasilkan bentuk gelombang tes dan menginterpretasikan hasil - Anda harus berhati-hati tentang itu - fase direferensikan ke waktu = 0, yaitu ketika Anda memulai pengambilan sampel (yaitu: ketika Anda mengumpulkan r [0] ), dan cosinus adalah interpretasi yang benar).

Kode di atas tidak elegan atau efisien (misalnya: gunakan tabel pencarian untuk nilai dosa / cos, dll.).

Hasil Anda akan menjadi lebih akurat karena Anda menggunakan N yang lebih besar, dan ada sedikit kesalahan karena fakta bahwa laju sampel dan N di atas tidak kelipatan satu sama lain.

Tentu saja, jika Anda ingin mengubah laju sampel, N, atau f, Anda harus mengubah kode dan nilai k. Anda dapat memasukkan nampan DFT ke mana saja di jalur frekuensi kontinu - pastikan Anda menggunakan nilai k yang sesuai dengan frekuensi yang diinginkan.

Masalahnya dapat dirumuskan sebagai masalah (nonlinier) kuadrat-terkecil:

$F(\phi)$$\phi$

Turunannya sangat sederhana:

$F''(\phi)$

$\phi^{2}$$2\pi k, k\in N$

Ada beberapa formulasi berbeda dari algoritma Goertzel. Yang menyediakan 2 variabel keadaan (ortogonal atau dekat dengan), atau variabel keadaan kompleks, karena keluaran yang mungkin sering dapat digunakan untuk menghitung atau memperkirakan fase dengan merujuk ke beberapa titik di jendela Goertzel, seperti tengah. Yang menyediakan output skalar tunggal saja biasanya tidak bisa.

Anda juga perlu tahu di mana jendela Goertzel Anda terkait dengan sumbu waktu Anda.

Jika sinyal Anda bukan bilangan bulat periodik di jendela Goertzel Anda, perkiraan fase di sekitar titik referensi di tengah jendela mungkin lebih akurat daripada merujuk fase ke awal atau akhir.

FFT penuh berlebihan jika Anda tahu frekuensi sinyal Anda. Plus Goertzel dapat disetel ke frekuensi yang tidak periodik dalam panjang FFT, sedangkan FFT akan membutuhkan interpolasi tambahan atau zero padding untuk frekuensi non-periodik-dalam-jendela.

Goertzel kompleks setara dengan 1 nampan DFT yang menggunakan pengulangan untuk vektor-vektor berbasis cosinus dan sinus atau faktor-faktor FFT.

Jika sinyal Anda bebas noise, Anda dapat mengidentifikasi nol penyilangan di keduanya dan menentukan frekuensi dan fase relatif.

$\phi$

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

Begitu:

• Apa yang Anda maksud dengan "puasa"?
• Seberapa akurat Anda membutuhkan perkiraan?
• $\phi$
• Berapa tingkat kebisingan pada setiap sinyal?

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

Untuk mengubah sinyal diskrit integral ke jumlah dan hati-hati memilih T!

Anda juga bisa melakukan ini (dalam notasi numpy):

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

di mana sinyal adalah sinyal fase-shifted Anda, cos dan sin adalah sinyal referensi, dan Anda menghasilkan perkiraan integral selama waktu tertentu melalui penjumlahan dari dua produk.

Ini merupakan peningkatan atas saran @Kevin McGee untuk menggunakan DFT frekuensi tunggal dengan indeks bin fraksional. Algoritma Kevin tidak membuahkan hasil yang luar biasa: sementara di setengah bins dan seluruh bins sangat tepat, juga dekat dengan keseluruhan dan membagi dua itu juga cukup bagus, tetapi jika tidak error bisa dalam 5%, yang mungkin tidak dapat diterima untuk sebagian besar tugas .

$N$$k$$N$

Kode di bawah ini dalam Swift, tetapi harus jelas secara intuitif:

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