Hitung dan interpretasikan frekuensi sesaat

Saya baru mengenal prinsip penghitungan frekuensi sesaat, dan muncul dengan banyak pertanyaan. Anda menemukan semuanya dalam daftar poin di akhir teks ini. Teksnya mungkin agak panjang, maafkan saya untuk itu, tapi saya benar-benar mencoba untuk menyelesaikan masalah itu sendiri.

Jadi saya tertarik pada frekuensi sesaat $f(t)$ dari sinyal bernilai nyata $x(t)$ . Perhitungan dilakukan dengan bantuan sinyal analitik $z(t) = x(t) + j y(t)$ , di mana $y(t)$ adalah transformasi Hilbert dari .$x(t)$

Untuk menghitung frekuensi sesaat dari sinyal analitik saya mengikuti makalah:$z(t)$

Perhitungan frekuensi sesaat dan bandwidth sesaat oleh Arthur E. Barns dari tahun 1992. Dalam makalah ini ia memperkenalkan beberapa metode untuk menghitung frekuensi sesaat. Saya menulis, semua formula yang ia usulkan (dan saya gunakan) sebentar lagi.

Untuk "belajar", saya bermain-main dengan sinyal yang sangat sederhana, dan dua sinyal yang lebih kompleks, di MATLAB, dan ingin mendapatkan frekuensi sesaat.

Fs = 1000;                                            % sampling-rate = 1kHz
t = 0:1/Fs:10-1/Fs;                                    % 10s 'Timevector'
chirp_signal = chirp(t,0,1,2);                         % 10s long chirp-signal, signal 1
added_sinusoid = chirp_signal + sin(2*pi*t*10);        % chirp + sin(10Hz), signal 2
modulated_sinusoid = chirp_signal .* sin(2*pi*t*10);   % chirp * sin(10Hz), signal 3

Plot dalam domain waktu dari ketiga sinyal tersebut terlihat sebagai berikut:

Plot dari semua frekuensi instan yang saya dapatkan setelah menerapkan semua metode dari kertas, adalah sebagai berikut:

Frekuensi sesaat dari sinyal kicauan murni: Frekuensi sesaat dari sinyal kicauan dengan tambahan sinusoid: Frekuensi sesaat dari sinyal kicauan termodulasi: Harap dicatat, bahwa dalam ketiga gambar, sumbu y plot 3 dan 4 diperbesar, sehingga amplitudo dari sinyal kicauan murni sinyal sangat kecil!

Kemungkinan pertama untuk mendapatkan sinyal analitik ke frekuensi sesaat adalah:

$$\begin{equation} f_{2}(t) = \frac{1}{2\pi} \frac{d}{dt}\theta(t) \end{equation}$$ $\theta(t)$

function [instantaneous_frequency] = f2(analytic_signal,Fs)
    factor =  Fs/(2*pi);
    instantaneous_frequency = factor * diff(unwrap(angle(analytic_signal)));
    % Insert leading 0 in return-vector to maintain size
    instantaneous_frequency = [0 instantaneous_frequency];
end

Dalam makalah, Lumbung sekarang menyarankan (atau lebih tepatnya mengatakan kompilasi) empat cara lain untuk menghitung frekuensi sesaat dari sinyal analitik. Dia juga menyebutkan formula atas, tetapi berpendapat bahwa itu tidak praktis karena ambiguitas dalam fase. Saya kira, dia tidak tahu unwrap()metodenya, atau lebih tepatnya matematika di belakangnya. (Saya sendiri baru tahu tentang metode itu hari ini, ketika melihat beberapa kode sumber lain pada frekuensi sesaat)

Dalam makalahnya, rumus memiliki label Nomor (2), oleh karena itu, saya memberikan f (t) indeks 2. Semua indeks lainnya sesuai dengan cara yang sama dengan nomor mereka di kertas.

Karena ambiguitas dalam fase, ia lebih suka menyarankan:

$$\begin{equation} f_{3}(t) = \frac{1}{2\pi}\cdot\frac{\overbrace{x(t)y'(t)}^\text{a}-\overbrace{x'(t)y(t)}^{\text{b}}}{\underbrace{x(t)^2}_{\text{c}}+\underbrace{y(t)^2}_{\text{d}}} \end{equation}$$

function [instantaneous_frequency] = f3(analytic_signal,Fs,T)
    x = real(analytic_signal);
    y = imag(analytic_signal);
    diff_x = diff(x);
    diff_y = diff(y);
    factor = Fs/(2*pi);
    a = x(2:end).*diff_y;
    b = y(2:end).*diff_x;
    c = x(2:end).^2;
    d = y(2:end).^2;
    instantaneous_frequency = factor * ((a-b)./(c+d));
    % Insert leading 0 in return-vector to maintain size
    instantaneous_frequency = [0 instantaneous_frequency];
end

$$\begin{equation} f_{9}(t) = \frac{1}{2\pi T}\cdot\text{arctan}\left[\frac{\overbrace{x(t)y(t+T)}^{\text{a}}-\overbrace{x(t+T)y(t)}^{\text{b}}}{\underbrace{x(t)x(t+T)}_{\text{c}}+\underbrace{y(t)y(t+T)}_{\text{d}}}\right] \end{equation}$$

function[instantaneous_frequency] = f9(analytic_signal, Fs, T)
    x = real(analytic_signal);
    y = imag(analytic_signal);
    factor = Fs/(2*pi*T);
    a = x(1:end-T).*y(1+T:end);
    b = x(1+T:end).*y(1:end-T);
    c = x(1:end-T).*x(1+T:end);
    d = y(1:end-T).*y(1+T:end);
    instantaneous_frequency = factor.*atan((a-b)./(c+d));
    % Append 0 to return-vector to maintain size
    instantaneous_frequency = [instantaneous_frequency zeros(1,T)];
end

$$\begin{equation} f_{11}(t) = \frac{1}{4\pi T}\cdot\text{arctan}\left[\frac{\overbrace{x(t-T)y(t+T)}^{\text{a}}-\overbrace{x(t+T)y(t-T)}^{\text{b}}}{\underbrace{x(t-T)x(t+T)}_{\text{c}}+\underbrace{y(t-T)y(t+T)}_{\text{d}}}\right] \end{equation}$$

function [instantaneous_frequency] = f11(analytic_signal, Fs, T)
    x = real(analytic_signal);
    y = imag(analytic_signal);
    factor = Fs/(4*pi*T);
    a = x(1:end-2*T).*y(1+2*T:end);
    b = x(1+2*T:end).*y(1:end-2*T);
    c = x(1:end-2*T).*x(1+2*T:end);
    d = y(1:end-2*T).*y(1+2*T:end);
    instantaneous_frequency = factor.*atan((a-b)./(c+d));
    % Append and insert 0s to maintain size
    instantaneous_frequency = [zeros(1,T) instantaneous_frequency zeros(1,T)];
end

$$\begin{equation} f_{14}(t)=\frac{2}{\pi T}\left[\frac{\overbrace{x(t)y(t+T)}^{\text{a}}-\overbrace{x(t+T)y(t)}^{\text{b}}}{\underbrace{(x(t)+x(t+T))^2}_{\text{c}}+\underbrace{(y(t)+y(t+T))^2}_{\text{d}}}\right] \end{equation}$$

function [instantaneous_frequency] = formula14(analytic_signal, Fs, T);
    x = real(analytic_signal);
    y = imag(analytic_signal);
    factor = 2*Fs/(pi*T);
    a = x(1:end-T).*y(1+T:end);
    b = x(1+T:end).*y(1:end-T);
    c = (x(1:end-T)+x(1+T:end)).^2;
    d = (y(1:end-T)+y(1+T:end)).^2;
    instantaneous_frequency = factor * ((a-b)./(c+d));
    % Append and insert 0s to maintain size
    instantaneous_frequency = [instantaneous_frequency zeros(1,T)];
end

Dalam semua 3 rumus perkiraan T diatur ke Fs (T = Fs = 1000 = 1s), seperti yang disarankan dalam makalah.

Sekarang pertanyaan saya adalah:

• Rumus f2 dan f3 mengembalikan hasil yang sama untuk sinyal kicauan murni. Saya pikir itu bagus, karena mereka menghitung sama. Tiga metode pendekatan tidak mengembalikan yang sama, bahkan sesuatu yang dekat dengannya! Kenapa begitu? (Saya harap ini bukan hanya bug pemrograman ...)
• Meskipun mereka kembali sama, terutama pada akhir plot mereka mulai 'menggoyangkan' banyak . Apa penjelasan untuk itu? Saya pertama kali memikirkan sesuatu seperti aliasing, tetapi frekuensi sampling saya cukup tinggi, dibandingkan dengan frekuensi sinyal, jadi saya pikir itu dapat dikecualikan.

• Setidaknya f2 dan f3 tampaknya bekerja sesuai pada sinyal kicauan murni, tetapi semua metode, termasuk f2 dan f3 tampaknya gagal mengerikan, ketika menyangkut lebih dari satu frekuensi dalam sinyal. Pada kenyataannya memiliki lebih dari satu frekuensi dalam suatu sinyal agak selalu demikian. Jadi bagaimana seseorang bisa mendapatkan (lebih atau kurang) frekuensi sesaat yang benar?

  • Sebenarnya saya bahkan tidak tahu apa yang diharapkan, ketika lebih dari satu frekuensi hadir dalam sinyal. Perhitungan mengembalikan satu angka untuk titik waktu tertentu, jadi apa yang harus dilakukan ketika, seperti di sini, lebih banyak frekuensi hadir? Mengembalikan rata-rata semua frekuensi atau sesuatu seperti itu?

• Dan pertanyaan saya yang mungkin paling penting adalah, bagaimana hal itu ditangani dalam perangkat lunak yang nyata dan rumit? Katakanlah saya ingin mengetahui frekuensi sesaat dari sinyal termodulasi pada 1,75 s, dan saya memilih metode f2, daripada saya bisa 'beruntung' dan mendapatkan angka mendekati 6 [Hz] yang kemungkinan besar merupakan jawaban yang benar, atau saya ambil hasil saya beberapa sampel di sebelahnya dan tiba-tiba saya mendapatkan beberapa kabel, cara ke tinggi, hasil, karena saya sayangnya memilih nilai di spike. Bagaimana ini bisa ditangani? Dengan postprocessing dengan rata-rata atau bahkan lebih baik filter median? Saya pikir bahkan itu mungkin menjadi sangat sulit terutama di daerah di mana banyak paku bersebelahan.

Dan satu pertanyaan terakhir, yang tidak begitu penting, mengapa sebagian besar makalah yang saya temukan pada frekuensi sesaat berasal dari area geografi, terutama dalam menghitung peristiwa seismografik seperti gempa bumi. Makalah Barne juga mengambil itu sebagai contoh. Bukankah frekuensi sesaat menarik di banyak daerah?

Sejauh ini, saya sangat berterima kasih atas setiap balasan, terutama ketika seseorang memberi saya tips tentang cara mengimplementasikannya dalam proyek perangkat lunak nyata ;)

Salam, Patrick

Tidak benar-benar jawaban tetapi mungkin membantu: Secara pribadi saya menemukan bahwa konsep frekuensi sesaat hanya berguna untuk sinyal band yang cukup sempit.

Pertimbangkan contoh sederhana dari dua gelombang sinus stabil, katakanlah 100Hz dan 934Hz. Dalam hal ini Anda tentu dapat menentukan dan menghitung frekuensi sesaat (dengan cara apa pun yang Anda inginkan) tetapi apa hasilnya? Wawasan atau properti apa yang mungkin dapat dimiliki oleh frekuensi sesaat yang mengatakan sesuatu yang bermakna tentang sinyal? Menerapkan konsep frekuensi sesaat ke sinyal yang memiliki banyak frekuensi pada saat yang sama, sama sekali tidak masuk akal.

Itu sebabnya Anda mendapatkan hasil yang layak untuk sapuan tetapi kurva aneh untuk sapuan + sinus. Ini juga merupakan alasan mengapa Anda melihat goyangan bagian tinggi dari sapuan. Bandwidth sinyal menjadi terlalu tinggi untuk menetapkan nomor frekuensi tunggal sehingga hasilnya melonjak.

Setidaknya f2 dan f3 tampaknya bekerja sesuai pada sinyal kicauan murni, tetapi semua metode, termasuk f2 dan f3 tampaknya gagal mengerikan, ketika menyangkut lebih dari satu frekuensi dalam sinyal. Pada kenyataannya memiliki lebih dari satu frekuensi dalam suatu sinyal agak selalu demikian. Jadi bagaimana seseorang bisa mendapatkan (lebih atau kurang) frekuensi sesaat yang benar?

seperti yang Hilmar sarankan, metode Hilbert transform (atau "Analytic Signal") tidak bekerja pada pita lebar karena ada lebih dari satu komponen frekuensi. Anda dapat melakukan metode ini hanya untuk komponen sinusoidal tunggal.

jadi, dengan pendekatan Sinyal Analitik, yang ingin Anda lakukan adalah memanfaatkan identitas ini:

$$\arctan u - \arctan v = \arctan \left( \frac{u-v}{1+uv} \right)$$

$|u-v|$$f_9$

tetapi harus ada hanya satu sinusoid yang berbeda-waktu dalam perhitungan transformasi Hilbert untuk melakukannya dengan benar. dan Anda lebih baik berbaris komponen "in-phase" dengan output dari transformasi Hilbert (yang ditunda dengan filter FIR kausal). kalau tidak, Anda akan mendapatkan omong kosong.

Wow, pertanyaan yang sangat besar. Saya akan menjawab pertanyaan yang tidak begitu penting terlebih dahulu:

Dan satu pertanyaan terakhir, yang tidak begitu penting, mengapa sebagian besar makalah yang saya temukan pada frekuensi sesaat berasal dari area geografi, terutama dalam menghitung peristiwa seismografik seperti gempa bumi. Makalah Barne juga mengambil itu sebagai contoh. Bukankah frekuensi sesaat menarik di banyak daerah?

Alasannya adalah bahwa sistem seismografik "vibroseis" digunakan dalam industri minyak untuk melakukan survei seismik. Truk yang saya tautkan bergetar dari sekitar 5 Hz hingga sekitar 90 Hz dan dapat dibuat untuk melakukan sinyal kicauan. Ada banyak uang di industri minyak, dan memproses pengembalian dari sinyal-sinyal ini bisa sangat, sangat menguntungkan. Oleh karena itu, banyak orang telah menghabiskan banyak waktu menganalisis sinyal-sinyal seperti itu, termasuk melihat teknik frekuensi sesaat.

---

$^{\rm TM}$

Lihatlah makalah ini.

Pendekatan yang lebih baik cenderung menggunakan "rata-rata tertimbang fase" seperti yang diterapkan di sini . Atau di sini untuk tautan langsung ke matlab .

Maaf memberikan jawaban setahun setelah fakta, tapi saya sengaja menemukan posting ini saat mencari artikel tentang topik ini. Pertanyaan Anda mencerminkan perselisihan dan interpretasi yang tersebar luas tentang "frekuensi sesaat" yang telah mengganggu bidang ini sejak awal. Banyak orang akan memberi tahu Anda, seperti beberapa jawaban di sini, bahwa IF hanya berlaku untuk sinyal "narrowband" atau "mono-komponen". Bahkan, itu tidak benar: kadang-kadang IF yang diperoleh oleh transformasi Hilbert berperilaku baik untuk sinyal broadband dan / atau "multi-komponen". Satu kuantitas yang telah diusulkan yang menghindari banyak kesulitan ini adalah "weighted average instantaneous frequency (WAIF)", yang dapat diukur menggunakan spektogram.

Lihat Loughlin di J. Acoust. Soc. Am., Januari 1999. Makalah bagus lainnya tentang IF dan kesalahpahaman umum adalah oleh Picinbono (IEEE Trans. Sig. Proc., Maret 1997) dan Vakman (IEEE Trans. Sig. Proc., April 1996).