👨🏻‍✈️ 🐓 🏂🏿 Transformasi Fourier. Cepat dan geram 🚴🏼 🤷 🧑🏽‍🤝‍🧑🏽

Seringkali, ketika mengembangkan algoritma, kita mengalami batas kompleksitas komputasi, yang, tampaknya, tidak mungkin diatasi. Transformasi Fourier memiliki kompleksitas

O (n^{2})

$O (n ^ 2)$ , dan versi cepat, diusulkan sekitar 1805 oleh House ¹ (dan diciptakan kembali pada tahun 1965 oleh James Cooley dan John Tukey)

O (n l o g (n))

$O (nlog (n))$ . Pada artikel ini saya ingin menunjukkan kepada Anda bahwa Anda bisa mendapatkan hasil konversi dalam waktu linier

O (n)

$O (n)$ atau bahkan mencapai kesulitan yang konstan

O (1)

$O (1)$ dalam kondisi tertentu yang ditemukan dalam masalah nyata.
Analisis Fourier

Ketika saya dihadapkan dengan tugas menulis program untuk menganalisis fungsi transfer sistem suara secara real time, saya, seperti orang lain, pertama-tama beralih ke konversi cepat. Semuanya baik-baik saja, tetapi dengan ukuran jendela waktu yang besar, beban CPU menjadi besar dan sesuatu harus dilakukan. Diputuskan untuk berhenti dan mempelajari transformasi lagi, dan pada saat yang sama mencari cara untuk menyelesaikan masalah. Mari kita kembali ke transformasi asli Joseph Fourier ² :

f (x) = j u m l a h b a t a s a n_{- i n f t y}^{+ i n f t y} c_{k} e^{2 p i i k x / T} c_{k} = f r a c 1 T i n t l i m i t_{0}^{T} f (x) e^{- 2 p i i k x / T} d x

$f (x) = \ jumlah \ batasan _ {- \ infty} ^ {+ \ infty} c_ke ^ {2 \ pi ikx / T} \\ c_k = \ frac {1} {T} \ int \ limit_0 ^ Tf ( x) e ^ {- 2 \ pi ikx / T} dx$

Kami akan melihat dengan cermat apa yang terjadi di sini. Setiap nilai output di domain frekuensi

c_{k}

$c_k$ adalah jumlah dari semua nilai input sinyal

f (x)

$f (x)$ dikalikan dengan

e^{- 2 p i i k x / T}

$e ^ {- 2 \ pi ikx / T}$ . Untuk membuat perhitungan, kita harus melalui semua data input untuk setiap nilai output, mis. untuk memenuhi itu

n^{2}

$n ^ 2$ operasi.

Singkirkan n

Biarkan saya mengingatkan Anda bahwa awalnya tugasnya adalah menganalisis data suara secara real time. Untuk melakukan ini, jendela waktu yang dipilih (dasarnya penyangga) ukuran N diisi dengan data dengan frekuensi yang sesuai dengan frekuensi pengambilan sampel. Dengan periode T, data input dikonversi dari jendela waktu ke jendela frekuensi. Jika Anda melihat bilangan real, maka N bervariasi dari 2 ¹⁴ (16 384) hingga 2 ¹⁶ (65 536) sampel (nilai diwarisi dari FFT, di mana ukuran jendela haruslah kekuatan dua). Waktu T = 80ms (12,5 kali per detik), yang memungkinkan Anda melihat perubahan dengan cukup nyaman dan tidak membebani CPU dan GPU. Frekuensi pengambilan sampel f _{d adalah} standar dan 48 kHz. Mari kita hitung berapa banyak data di jendela waktu yang berubah antar dimensi. Selama waktu T, ia memasuki buffer

48000 * 0, 08 = $ 384

$48000 * 0,08 = $ 384$ sampel. Dengan demikian, hanya 5% hingga 23% dari data diperbarui di jendela. Dalam kasus terburuk, 95% (dan paling banyak 73%, yang juga banyak!) Dari sampel yang diproses akan jatuh ke dalam konversi lagi dan lagi, meskipun fakta bahwa mereka telah diproses dalam iterasi sebelumnya.

Pembaca yang penuh perhatian pada saat itu akan mengangkat tangannya dan berkata: "tunggu, tapi bagaimana dengan koefisien

e^{- 2 p i i k x / T}

$e ^ {- 2 \ pi ikx / T}$ ? Lagi pula, dengan setiap transformasi baru, data yang sama akan ditempatkan di posisi baru seri, dan sebagai hasilnya, memiliki koefisien yang berbeda? " Untuk setiap lima untuk perawatan mereka, mari kita ingat detail penting dalam transformasi yang sering dilupakan. Dalam studi tentang nilai-nilai fungsi

f (t)

$f (t)$ selama interval dari 0 hingga t, fungsi tersebut dianggap periodik, yang memungkinkan Anda untuk menggeser fungsi ke kiri atau kanan dalam waktu tanpa kesulitan. Akibatnya, kami berhak untuk tidak memasukkan nilai baru di akhir dan menghapus nilai lama dari awal, tetapi untuk mengganti data secara siklikal dalam buffer.

Untuk kejelasan, Anda bisa menulis dalam bentuk tabel bagaimana buffer akan berubah:

t = 0	f (0)	f (1)	f (2)	f (3)	f (4)	f (5)	f (6)	f (7)	f (8)	f (9)
t = 1	f (10)	f (1)	f (2)	f (3)	f (4)	f (5)	f (6)	f (7)	f (8)	f (9)
t = 2	f (10)	f (11)	f (2)	f (3)	f (4)	f (5)	f (6)	f (7)	f (8)	f (9)
t = 3	f (10)	f (11)	f (12)	f (3)	f (4)	f (5)	f (6)	f (7)	f (8)	f (9)
t = 4	f (10)	f (11)	f (12)	f (13)	f (4)	f (5)	f (6)	f (7)	f (8)	f (9)

Anda dapat menulis bagaimana transformasi dalam waktu berubah dari t ₁ ke t ₂ :

F_{t} = F_{t - 1} + D e l t a F D e l t a F : D e l t a c_{k} = f r a c 1 T i n t l i m i t_{t_{1}}^{t_{2}} (f_{t} (x) - f_{t - 1} (x)) e^{- 2 p i i k x / T} d x

$F_t = F_ {t-1} + \ Delta F \\ \ Delta F: \ Delta c_k = \ frac {1} {T} \ int \ limit_ {t_1} ^ {t_2} (f_t (x) - f_ { t-1} (x)) e ^ {- 2 \ pi ikx / T} dx$

Nilai

F_{t - 1} (x)

$F_ {t-1} (x)$ adalah hasil dari konversi sebelumnya, dan kompleksitas perhitungan

D e l t a f (x)

$\ Delta f (x)$ tidak tergantung pada ukuran jendela waktu dan, oleh karena itu, konstan. Akibatnya, kompleksitas konversi akan menjadi

O (n)

$O (n)$ ^* karena semua yang tersisa bagi kita adalah melalui jendela frekuensi sekali dan menerapkan perubahan untuk sampel T yang telah berubah selama perjalanan waktu. Saya juga ingin menarik perhatian Anda pada fakta bahwa peluangnya

e^{- 2 p i i k x / T}

$e ^ {- 2 \ pi ikx / T}$ dapat dihitung terlebih dahulu, yang memberikan tambahan produktivitas, dan hanya dua operasi yang tersisa dalam siklus: mengurangi bilangan real dan mengalikan bilangan real dengan yang kompleks, dalam praktiknya kedua operasi ini sederhana dan murah.

Untuk melengkapi gambar, tetap hanya untuk menunjukkan keadaan awal, tetapi di sini semuanya sederhana:

F_{0} (x) = 0

$F_0 (x) = 0$

* - tentu saja, kompleksitas akhir dari seluruh transformasi akan tetap demikian

O (n^{2})

$O (n ^ 2)$ , tetapi akan dieksekusi secara bertahap, lebih dari n iterasi, sementara buffer diperbarui.

O (n)

$O (n)$ - ini adalah kompleksitas memperbarui data, tetapi ini persis apa yang kita butuhkan (saat menggunakan FFT, kompleksitas setiap transformasi

O (n l o g (n))

$O (nlog (n))$ )

Tetapi bagaimana jika Anda menggali lebih dalam. Atau singkirkan yang kedua n

Saya ingin membuat reservasi segera bahwa langkah-langkah selanjutnya hanya berlaku jika Anda tidak berencana untuk melakukan transformasi terbalik untuk hasilnya (untuk memperbaiki sinyal atau mendapatkan tanggapan impuls). Untuk memulainya, saya ingin mengingatkan Anda bahwa sebagai hasil konversi, kami mendapatkan larik data simetris, yang segera memungkinkan kami untuk mengurangi jumlah konversi hingga setengahnya.

Sekarang mari kita menganalisis set data yang dihasilkan, mengingat kondisi masalah. Kami memiliki satu set bilangan kompleks, yang masing-masing menggambarkan amplitudo dan fase osilasi pada frekuensi tertentu. Frekuensi dapat ditentukan dengan rumus:

f [j] = j f r a c f d N

$f [j] = j \ frac {fd} {N}$ untuk

j < f r a c N 2

$j <\ frac {N} {2}$ . Mari mengevaluasi langkah jendela frekuensi pada data kami:

D e l t a f = f r a c f d N

$\ Delta f = \ frac {fd} {N}$ Untuk N = 2 ¹⁴ : 2.93 Hz (dan untuk 2 ¹⁶ : 0.73 Hz). Dengan demikian, dalam kisaran dari 1 kHz hingga 2 kHz kami mendapatkan 341 hasil. Cobalah untuk mengevaluasi secara independen berapa banyak data akan berada dalam kisaran dari 8 kHz hingga 16 kHz untuk N = 65536. Banyak, kan? Banyak! Apakah kita membutuhkan begitu banyak data? Tentu saja, dalam masalah menampilkan karakteristik frekuensi sistem suara, jawabannya adalah tidak. Dan di sisi lain, untuk analisis di wilayah frekuensi rendah, langkah kecil sangat berguna. Jangan lupa bahwa masih ada jadwal di depan yang harus volume ini (

f r a c N 2

$\ frac {N} {2}$ ) mengonversi ke bentuk yang dapat dibaca manusia (rata-rata, spline, atau menghaluskan) dan menampilkannya di layar. Dan pada frekuensi tinggi, bahkan memiliki layar 4K dan menampilkan grafik dalam mode layar penuh dengan sumbu frekuensi logaritmik, ukuran langkah dengan cepat berubah menjadi kurang dari 1 piksel.

Berdasarkan pengalaman, Anda dapat mengetahui bahwa cukup memiliki hanya 48 poin per oktaf, dan memiliki data sedikit lebih halus dan lebih rata-rata, saya sarankan berhenti di 96. Dalam rentang frekuensi audio dari 20 Hz hingga 20 kHz, mudah untuk menghitung hanya 10 oktaf: 20, 40, 80 , 160, 320, 640, 1280, 2560, 5120, 10240, 20480, yang masing-masing dapat dibagi menjadi sejumlah subband tertentu (jangan lupa bahwa partisi harus dilakukan secara geometris, dan bukan secara hitung), oleh karena itu, lebih dari cukup untuk melakukan konversi hanya untuk 960 frekuensi untuk didapatkan Performan yang di 16 ... 65 kali lebih kecil dari versi aslinya.

Dengan demikian, menggabungkan kedua pendekatan, kami memperoleh kompleksitas konstan dari algoritma pembaruan data

O (1)

$O (1)$ .

Madu kuadrat dan sesendok tar

Sekarang kita dapat mengatakannya dengan aman dari kerumitan

O (n^{2})

$O (n ^ 2)$ kami sampai pada kompleksitas

O (1)

$O (1)$ menggunakan dua peretas sederhana:

Setelah menganalisis masalah, kami perhatikan bahwa data ditambahkan secara bertahap, dan periode pembaruan lengkap dari time time jauh lebih tinggi daripada periode transformasi dan selanjutnya menghitung perbedaan dari transformasi Fourier.
beralih dari langkah aritmatika di jendela frekuensi menjadi hanya dibatasi oleh nilai yang ditentukan, yang secara dramatis dapat mengurangi jumlah konversi.

Tapi, tentu saja, hidup benar-benar akan menjadi dongeng, jika bukan satu melainkan. Penerapan kedua pendekatan ini memungkinkan kami untuk benar-benar membongkar CPU sehingga dapat menebak bahwa ia menghitung transformasi Fourier dan menampilkan hasil di layar bahkan dengan

N = 2^{20}

$N = 2 ^ {20}$ itu sulit. Tetapi hukuman tidak membuat dirinya menunggu ketika sinyal Anda pada kenyataannya tidak periodik (dan ini diperlukan untuk mendapatkan hasil konversi yang benar) dan tidak mungkin untuk memilih ukuran jendela yang sesuai, menjadi perlu untuk menggunakan berbagai fungsi jendela, yang tidak lagi memungkinkan Anda untuk menggunakan langkah pertama sepenuhnya. Praktek telah menunjukkan bahwa penggunaan fungsi jendela sangat penting dalam studi sinyal dengan frekuensi kurang dari

0, 1 f_{d}

$0,1f_d$ . Pada frekuensi tinggi, jumlah periode yang jatuh ke jendela waktu secara signifikan mengurangi distorsi yang timbul dari adanya celah tingkat pertama (antara f (0) dan f (N-1)) dalam fungsi asli.

Pada akhirnya, saya juga menolak langkah kedua dan kembali ke FFT, karena keuntungan dalam tugas ini sudah kecil.

Kesimpulannya

Pendekatan pertama dapat diterapkan jika data Anda bersifat periodik yang diucapkan dan perlu dianalisis dari waktu ke waktu menggunakan jendela waktu besar, yang, saya ingat, tidak harus berupa derajat 2, yaitu. nomor alami.
Pendekatan kedua dapat diterapkan (bahkan dengan mempertimbangkan fungsi jendela akun) jika hanya sejumlah kecil frekuensi yang dianalisis dalam data.

Sayangnya, bagi saya dalam masalah ini, ini hanya sedikit hiburan matematis, tapi saya harap itu menginspirasi Anda untuk mempelajari algoritma lain pada hari libur dalam hal perubahan dalam input data dari waktu ke waktu :)

Sastra

Gambar diambil dari manga oleh Michio Shibuya. “MATEMATIKA YANG MENARIK. Analisis Fourier

Transformasi Fourier. Cepat dan geram

Singkirkan n

Tetapi bagaimana jika Anda menggali lebih dalam. Atau singkirkan yang kedua n

Madu kuadrat dan sesendok tar

Kesimpulannya

Sastra

More articles: