Fungsi jendela DIY

Dalam pemrosesan sinyal digital, fungsi jendela banyak digunakan untuk membatasi sinyal dalam waktu, dan namanya dikenal baik oleh semua orang yang telah menemukan transformasi Fourier diskrit dalam satu atau lain cara: Hann, Hamming, Blackman, Harris, dan lainnya. Tetapi apakah mereka cukup, apakah mungkin untuk menghasilkan sesuatu yang baru dan apakah ada gunanya?



Pada artikel ini, kita akan melihat output dari fungsi jendela dengan properti baru menggunakan Wolfram Mathematica. Diasumsikan juga bahwa pembaca memiliki pemahaman umum tentang pemrosesan sinyal digital dalam konteks masalah yang sedang dibahas dan setidaknya akrab dengan artikel dari Wikipedia .

pengantar

Secara historis, fungsi jendela pertama muncul dalam proses upaya untuk meningkatkan sifat spektralnya dengan mengurangi amplitudo lobus samping - karena ketika sinyal dikalikan dengan fungsi jendela, spektrumnya berbelit-belit, yang menyebabkan beberapa batasan untuk analisis. Pada saat itu, tahun 50-an abad ke-20, kekuatan komputasi tidak memungkinkan untuk memanipulasi komputasi simbolik dengan mudah dan alami, dan pemilihan parameter yang optimal cukup sulit. Inilah salah satu alasan mengapa fungsi disebut dengan nama yang berbeda - mereka, atau lebih tepatnya, sifat spektralnya, telah dipelajari oleh peneliti yang berbeda untuk waktu yang cukup lama. Efek samping dari hal ini adalah bahwa kumpulan fungsi jendela bernama yang ada dianggap sebagai semacam "kumpulan standar" yang di luarnya tidak ada yang dapat ditemukan;pada saat yang sama, nama-nama jendela ini tidak membawa informasi apa pun tentang properti dan menyarankan studi dan penghafalan terpisah.

Klasifikasi

Semua fungsi jendela diperoleh dengan mengalikan beberapa fungsi dengan fungsi persegi panjang, yang dijamin akan memastikan batasan waktunya. Oleh karena itu, pertama-tama, mereka dapat diklasifikasikan menurut fungsi yang mereka gunakan:

1. jumlah cosinus. Kelas paling luas karena spektrumnya mudah dihitung dan merupakan jumlah dari fungsi sinc berbobot. Ini termasuk Hann, Hamming, Blackman, Harris ...
2. polinomial kontinu sedikit demi sedikit. Mereka diperoleh sebagai hasil dari konvolusi fungsi sederhana - misalnya, persegi panjang dan segitiga. Pada saat yang sama, spektrum mereka berlipat ganda dan penemuan mereka juga tidak terlalu sulit. Ini termasuk Dirichlet, Triangular, Parzen, Welch ...
3. semua yang lain - menggunakan eksponensial, Gaussians, sinc dan lain-lain, pilihan fungsi tertentu yang lebih bersifat ideologis daripada sifat spektral secara spesifik.

Mereka juga dapat dibagi berdasarkan properti di tepinya:

1. tidak ada celah di tepinya - nilai sama dengan nol. Istirahat adalah Dirichlet, Hamming, Blackman-Harris. Tidak punya - Segitiga, Hann, Nuttal
2. tidak adanya diskontinuitas dari turunan pertama dan lainnya

Dua fungsi jendela lainnya menonjol secara terpisah:

1. Kaiser, yang memungkinkan Anda mengatur lebar kelopak utama;
2. Dolph-Chebyshev, semua lobus sampingnya sama dengan amplitudo yang diberikan.

Keduanya memiliki diskontinuitas di tepi kedua nilai dan turunannya, dan perhitungannya melibatkan beberapa kesulitan - fungsi Kaiser dihitung melalui fungsi khusus (Bessel), dan fungsi Dolph-Chebyshev ditentukan dalam domain frekuensi dan dihitung melalui transformasi Fourier diskrit terbalik. Menemukan antiturunan analitis mereka juga sangat sulit.

Rekayasa terbalik

Mari kita lihat definisi dari fungsi Blackman dan Nuttel:

BlackmanWindow[x] // FunctionExpand

↓

$$

$$\left\{ \begin{array}{ll} \frac{1}{50} (25 \cos (2 \pi x)+4 \cos (4 \pi x)+21) & -\frac{1}{2}\leq x\leq \frac{1}{2} \\ 0 & \left| x\right| >\frac{1}{2} \\ \end{array} \right.$$

NuttallWindow[x] // FunctionExpand

↓

$$

$$\left\{ \begin{array}{ll} \frac{121849 \cos (2 \pi x)+36058 \cos (4 \pi x)+3151 \cos (6 \pi x)+88942}{250000} & -\frac{1}{2}\leq x\leq \frac{1}{2} \\ 0 & \left| x\right| >\frac{1}{2} \\ \end{array} \right.$$

Mereka adalah jumlah dari gelombang kosinus 3 dan 4 dengan frekuensi genap (mulai dari nol) dan beberapa koefisien. Dari mana koefisien ini berasal? Tidak mungkin mereka dipilih dengan tangan, setidaknya masing-masing secara terpisah - lagipula, ada kondisi batas yang harus dipatuhi jendela terlepas dari sifat spektralnya - setidaknya sama dengan kesatuan di tengah.



Mari kita coba "menemukan" fungsi Blackman sendiri. Untuk melakukan ini, kita mendefinisikan fungsi dengan koefisien yang belum diketahui

f = Function[x, a + b Cos[2 Pi x] + c Cos[4 Pi x]];

dan menyusun sistem persamaan yang mendefinisikan kondisi batas - persamaan ke kesatuan di tengah, nol di tepi, dan persamaan ke nol turunan di tepinya:

ss = Solve[{
   f[0] == 1,
   f[1/2] == 0,
   f'[1/2] == 0
   }, {a, b, c}]

↓

$$

$$\left\{\left\{b\to \frac{1}{2},c\to \frac{1}{2}-a\right\}\right\}$$

Solusinya ditemukan, tetapi tidak hanya satu, tetapi banyak - tergantung pada koefisien a , yang dengan sopan diperingatkan Wolfram kepada kita. Sekarang, dari solusi yang ditemukan, kita mendefinisikan fungsi baru tergantung pada koefisien yang tidak diketahui: ↓ Sangat mudah untuk melihat bahwa untuk a = 0,42 kita mendapatkan fungsi Blackman. Tapi kenapa tepatnya 0,42? Untuk menjawab pertanyaan ini, kita perlu memplot spektrumnya. Transformasi Fourier Analitik bukanlah tugas yang paling mudah, tetapi Wolfram dapat menanganinya juga, menyelamatkan kita dari tugas.



hx = Function[{x, a},

Piecewise[{{f[x] /. ss[[1]], -1/2 < x < 1/2}}, 0] // Evaluate]

hw = Function[{w, a}, 
  FourierTransform[hx[x, a], x, w] // #/Limit[#, w -> 0] & // 
    Simplify // Evaluate]

↓









Dalam bentuk ini, grafik dari fungsi tersebut masih belum terlalu informatif, karena lebih mudah untuk menganalisis spektrum pada skala logaritmik. Dengan menambahkan konversi yang sesuai ke desibel$$$20 \log _{10}(\left| x\right| )$ , kita mendapatkan grafik yang sama dari spektrum fungsi jendela:









Dalam proses mengubah parameter, kita akan mengamati hal yang serupa:







Tergantung pada parameter a, level lobus samping berubah, dan pada a = 0,42 itu lebih atau kurang minimal dan seragam. Dengan a = 0.409, kita bisa mendapatkan hasil yang sedikit lebih baik jika yang dimaksud dengan "lebih baik" adalah level minimal lobus samping.

Pengembangan

Jelas, pengoptimalan mikro seperti itu sama sekali tidak sepadan dengan usaha. Apakah mungkin untuk meningkatkan properti jendela ini secara kualitatif tanpa mengubah data awal?



Sebelumnya kita telah melihat output dari fungsi jendela yang berjumlah satu, yang memungkinkan Anda untuk mengembalikan sinyal asli. Mari kita coba teknik ini untuk jendela kita dan lihat bagaimana pengaruhnya terhadap spektrum.



Mari tentukan fungsi tambahan untuk membangun jendela yang tumpang tindih, dengan mempertimbangkan batasan argumen yang berubah dalam rentang dari -1/2 ke 1/2:







Temukan antiturunannya, geser ke pusat koordinat dan skala ke satu:





↓





tampilkan pada grafik:







Seperti yang Anda lihat, Wolfram melakukannya sendiri di sini juga, dan kami tidak harus secara manual menetapkan definisi antiturunan secara bertahap. Sekarang tampilan jendela kita tidak hanya bergantung pada variabel a , tetapi juga pada tingkat tumpang tindih - dan seiring bertambahnya, itu akan cenderung ke bentuk turunan:







Dan sentuhan terakhir adalah menemukan fungsi analitik untuk spektrum untuk menentukan nilai optimal dari parameter a .



Di sini, jika kita mencoba menghitung transformasi secara langsung, seperti terakhir kali, kita akan mengarahkan Wolfram ke pemikiran yang dalam. Ada beberapa cara untuk mempercepat proses ini:



- gunakan transformasi diskrit daripada transformasi analitis - yang paling sederhana. Tetapi ahli matematika sejati tidak akan menerapkan metode numerik di mana solusi analitis murni dapat ditemukan - dan kami tidak akan (belum); selain itu, ia memiliki batasan yang jelas (dalam hal resolusi dan batasan spektrum);



- berikan nilai tertentu sebagai parameter tumpang tindih dan, lihat hasilnya, coba lihat polanya dan generalisasikan. Pilihan yang cukup berhasil, tetapi membutuhkan upaya mental dan kreatif tambahan. Mari kita tinggalkan sebagai pilihan terakhir.



- membuat kalkulasi langsung dalam domain frekuensi. Hal ini dimungkinkan karena properti transformasi Fourier berikut:



1) linier, mis. jumlah gambar sama dengan gambar jumlah;

2) integrasi dalam domain waktu setara dengan membagi dengan iw dalam domain frekuensi.

3) perpindahan waktu setara dengan modulasi dengan sinusoid kompleks.

(lebih lanjut di wikipedia atau di sini ).



Jadi, dengan memiliki spektrum hw dari fungsi jendela sembarang, kita dapat memperoleh dari itu spektrum hwo dari suatu fungsi yang dijumlahkan menjadi satu dengan tumpang tindih t :





↓





Sekarang kita dapat melihat bagaimana spektrum berubah dalam dinamika:







Di sini, mengubah parameter tumpang tindih sudah akan mempengaruhi spektrum jendela dengan cara yang sedikit berbeda :







Setelah bermain sedikit dengan parameter, mudah untuk melihat bahwa "lebih banyak tidak lebih baik", dan tingkat optimal tumpang tindih untuk fungsi itu berada di wilayah empat. Khususnya untuk t = 4 dan a= 0,404 kita mendapatkan level lobus samping tidak melebihi -80 dB. Ini adalah hasil yang sangat bagus - terutama mengingat fungsi kosinus yang ditinggikan, yang secara tradisional digunakan untuk jendela penjumlahan, memberikan tingkat lobus sekitar -30 dB. Nah, ketika ditulis secara eksplisit, fungsi jendela baru kita akan terlihat seperti ini:

$$

$$\left\{ \begin{array}{ll} \frac{404 \pi (1-2 x)+36 \sin \left(\frac{1}{3} (16 \pi x+\pi )\right)+375 \sin \left(\frac{1}{3} (\pi -8 \pi x)\right)}{606 \pi } & \frac{1}{4} < x \leq \frac{1}{2} \\ \frac{404 (2 \pi x+\pi )+36 \sin \left(\frac{2}{3} \pi (8 x+1)\right)+375 \sin \left(\frac{1}{3} (8 \pi x+\pi )\right)}{606 \pi } & -\frac{1}{2} < x \leq -\frac{1}{4} \\ \frac{\sqrt{3} \left(125 \cos \left(\frac{8 \pi x}{3}\right)+12 \cos \left(\frac{16 \pi x}{3}\right)\right)}{202 \pi }+\frac{1}{3} & -\frac{1}{4} < x \leq \frac{1}{4} \\ 0 & \left| x\right| >\frac{1}{2} \\ \end{array} \right.$$

Pengembangan lebih lanjut

Apa lagi yang dapat Anda lakukan untuk mengurangi level sidelobe? Anda dapat mengambil cosinus bukan dengan genap, tetapi dengan frekuensi ganjil (di sini, untuk kekompakan, solusi sistem persamaan ditanamkan langsung ke dalam definisi fungsi):

$$

$$\left(a \cos (\pi x)+\left(\frac{5}{8}-\frac{a}{2}\right) \cos (3 \pi x)+\left(\frac{3}{8}-\frac{a}{2}\right) \cos (5 \pi x)\right) \Pi (x)$$

dan setelah integrasi dengan parameter a = 0,6628 dan tingkat tumpang tindih 4,5, dapatkan atenuasi -90 dB:







Secara eksplisit:

$$

$$\small\left\{ \begin{array}{ll} \frac{327 \sin \left(\frac{1}{7} \pi (45 x+2)\right)+24855 \sin \left(\frac{1}{7} \pi (1-9 x)\right)+3670 \cos \left(\frac{1}{14} \pi (54 x+1)\right)+21512}{43024} & \frac{5}{18} < x \leq \frac{1}{2} \\ \frac{24855 \sin \left(\frac{1}{7} \pi (9 x+1)\right)+3670 \sin \left(\frac{3}{7} \pi (9 x+1)\right)+327 \sin \left(\frac{5}{7} \pi (9 x+1)\right)+21512}{43024} & -\frac{1}{2} < x \leq -\frac{5}{18} \\ \begin{array}{ll} \frac{-3670 \sin \left(\frac{3}{7} \pi (9 x-1)\right)+24855 \sin \left(\frac{1}{7} \pi (9 x+1)\right)+327 \sin \left(\frac{5}{7} \pi (9 x+1)\right)}{43024}+ \\ +\frac{3670 \sin \left(\frac{3}{7} \pi (9 x+1)\right)+327 \sin \left(\frac{1}{7} \pi (45 x+2)\right)+24855 \sin \left(\frac{1}{7} \pi (1-9 x)\right)}{43024} \end{array} & -\frac{5}{18} < x \leq \frac{5}{18} \\ 0 & \left| x\right| >\frac{1}{2} \\ \end{array} \right.$$

Anda dapat menambahkan satu cosinus lagi dan menambah jumlah turunan nol:

$$

$$\left(a \cos (\pi x)+\frac{1}{80} (35-16 a) \cos (3 \pi x)+\frac{1}{80} (35-48 a) \cos (5 \pi x)+\frac{1}{40} (5-8 a) \cos (7 \pi x)\right) \Pi (x)$$

dan setelah integrasi dengan parameter a = 0,5862 dan tingkat tumpang tindih 6,4, dapatkan penekanan -110 dB:







Secara eksplisit:

$$

$$\small\left\{ \begin{array}{ll} \frac{-3077550 \sin \left(\frac{1}{54} \pi (64 x-5)\right)-90069 \sin \left(\frac{5}{54} \pi (64 x-5)\right)-5820 \sin \left(\frac{7}{54} \pi (64 x-5)\right)-560455 \cos \left(\frac{1}{9} \pi (7-32 x)\right)+2601344}{5202688} & \frac{11}{32}< x \leq \frac{1}{2} \\ \frac{560455 \sin \left(\frac{1}{18} \pi (64 x+5)\right)+3077550 \sin \left(\frac{1}{54} \pi (64 x+5)\right)+90069 \sin \left(\frac{5}{54} \pi (64 x+5)\right)+5820 \sin \left(\frac{7}{54} \pi (64 x+5)\right)+2601344}{5202688} & -\frac{1}{2}< x \leq -\frac{11}{32} \\ \begin{array}{ll} \frac{-3077550 \sin \left(\frac{1}{54} \pi (64 x-5)\right)-90069 \sin \left(\frac{5}{54} \pi (64 x-5)\right)-5820 \sin \left(\frac{7}{54} \pi (64 x-5)\right)-560455 \cos \left(\frac{1}{9} \pi (7-32 x)\right)}{5202688}+ \\ +\frac{560455 \sin \left(\frac{1}{18} \pi (64 x+5)\right)+3077550 \sin \left(\frac{1}{54} \pi (64 x+5)\right)+90069 \sin \left(\frac{5}{54} \pi (64 x+5)\right)+5820 \sin \left(\frac{7}{54} \pi (64 x+5)\right)}{5202688} \end{array} &-\frac{11}{32}< x \leq \frac{11}{32} \\ 0 & \left| x\right| >\frac{1}{2} \\ \end{array} \right.$$

Penurunan yang lebih signifikan dalam tingkat lobus samping dapat dicapai dengan mengkuadratkan transformasi Fourier dan dengan demikian mengurangi tingkat lobus samping dengan faktor 2 sekaligus. Hal ini memungkinkan Anda untuk menyingkirkan peningkatan jumlah parameter untuk pemilihan manualnya, tetapi menambah kerumitan dalam menghitung konvolusi dalam domain waktu.









dan di sini sudah mungkin untuk mencapai penekanan lebih dari 160 dB:







Kami tidak akan memberikan rumus secara eksplisit karena ukurannya yang mengesankan.

Fungsi jendela hipergeometrik

Untuk memastikan jumlah nol yang diperlukan di batas fungsi kami, kami menggunakan pencarian melalui solusi sistem persamaan dengan integrasi berikutnya. Ini sangat tidak nyaman - Anda harus mengubah jumlah persamaan setiap waktu. Mungkin ada solusi yang lebih sederhana dan lebih indah? Ada! Dan fungsi hipergeometrik akan membantu kita dalam hal ini.





Solusi kami akan terlihat seperti ini:

$$

$$\frac{2 z \Gamma \left(\frac{n+3}{2}\right) \, _2F_1\left(\frac{1}{2},-\frac{n}{2};\frac{3}{2};z^2\right)}{\sqrt{\pi } \Gamma \left(\frac{n}{2}+1\right)}+a \left(z+b z^3\right) \left(1-z^2\right)^{n/2}$$

Dimana $$$z=\sin \left(\frac{\pi x}{2}\right)$



Ini terdiri dari jumlah dua bagian: antiturunan fungsi $$$\left(1-z^2\right)^{\frac{n}{2}}$dinormalisasi dalam amplitudo ke titik (1,1) dan fungsi dengan parameter bebas, dikalikan dengan yang sama $$$\left(1-z^2\right)^{\frac{n}{2}}$untuk mempertahankan jumlah turunan nol yang diperlukan:







Keakuratan penyesuaian dan pengaruhnya terhadap distribusi lobus samping akan bergantung pada fungsi mana yang kita pilih untuk mengoreksi bentuk - ada kemungkinan opsi yang memerlukan studi terpisah. Dalam hal ini, fungsinya$$$a \left(z+b z^3\right)$memberikan opsi yang sepenuhnya dapat diterima - dan karena dua parameter memungkinkan Anda membuat pengaturan yang lebih tepat.



Rumusnya sendiri menarik (dan nyaman) karena untuk genap n disederhanakan menjadi polinomial, dan untuk ganjil n ke jumlah polinomial dengan busur dan akar kuadrat - yang membuat rumus akhir lebih kompak dan lebih mudah untuk dihitung.



Lebih mudah (dan lebih cepat) untuk memilih parameter di sini melalui transformasi Fourier diskrit. Kita juga membutuhkan definisi tambahan dari fungsi overlap untuk bekerja dengan tiga parameter:



















Sebagai contoh, setelah mengganti parameter dari gambar, fungsi kita akan disederhanakan menjadi

$$

$$\frac{-288 z^{11}+2315 z^9-7380 z^7+12330 z^5-11940 z^3+8163 z}{3200}$$

Jendela Kaiser Terbalik

Jika tugas menjumlahkan jendela menjadi satu tidak sepadan, maka jendela paling optimal adalah jendela Kaiser, yang memungkinkan Anda menyesuaikan lebar lobus utama dengan mulus. Namun, ini memiliki kelemahan - karena diekspresikan dalam fungsi Bessel , agak sulit untuk membacanya di luar paket matematika. Anda dapat, tentu saja, menerapkan fungsi Bessel secara terpisah - atau Anda dapat mencari perkiraan melalui fungsi dasar. Dan tiba-tiba ternyata dengan menggunakan fungsi spektrum jendela Kaiser (waktu terbatas) untuk ini, Anda bahkan bisa mendapatkan hasil yang sedikit lebih baik -

$$

$$\frac{ \sinh \left(k \sqrt{1-4 x^2}\right)}{\sinh (k) \sqrt{1-4 x^2} } \Pi (x)$$

Pembayaran untuk keputusan licik seperti itu adalah ketidakmungkinan mengekspresikan spektrumnya secara analitis - tetapi ini tidak begitu penting, karena dalam praktiknya, dalam kasus apa pun, kami mengoperasikan dan menganalisis data dalam bentuk diskrit menggunakan transformasi Fourier diskrit. Di bawah pada grafik Anda dapat membandingkan spektrumnya, di mana merah adalah jendela Kaiser asli, hijau adalah perkiraannya:









Dengan lebar kelopak utama yang sama, tingkat lobus samping sedikit lebih rendah. Dan untuk kemudahan penggunaan, Anda dapat membuat tabel dengan nilai perkiraan parameter k untuk mendapatkan tingkat penekanan sidelobe yang diperlukan:

$$

$$\begin{array}{cc} \text{-60 } & \text{k=8.8} \\ \text{-90 } & \text{k=11.36} \\ \text{-120 } & \text{k=15.18} \\ \text{-150 } & \text{k=18.88} \\ \end{array}$$

Menemukan perkiraan antiturunan dari jendela Kaiser terbalik tetap menjadi masalah menarik yang terpisah. Sejauh ini, hanya mungkin untuk mengekspresikannya melalui deret pangkat, yang menyatu agak lambat:

$$

$$\small \int \frac{ \sinh \left(k \sqrt{1-4 x^2}\right)}{\sinh (k) \sqrt{1-4 x^2} } \, dx= \sum _{n=1}^{\infty } -\frac{2^{\frac{1}{2}-n} k^{n-1} x^{2 n-1} \left(\sqrt{-k} (-1)^n K_{n-\frac{1}{2}}(-k)+\sqrt{k} K_{n-\frac{1}{2}}(k)\right)}{\sqrt{\pi } (2 n-1) \sinh (k) \Gamma (n)}$$

Mungkin para ahli matematika dapat menyarankan solusi yang lebih baik di kolom komentar.

Kesimpulan

Terlepas dari kenyataan bahwa pemrosesan sinyal digital adalah disiplin yang sepenuhnya terbentuk, masih ada ruang untuk pengembangan dan menemukan sesuatu yang baru di dalamnya - setidaknya dalam literatur khusus dan ilmiah yang ada, masalah membangun fungsi jendela yang diringkas menjadi satu unit tidak dipertimbangkan; dan kapabilitas sistem aljabar komputer modern memungkinkan pemikiran ulang pengetahuan yang terkumpul pada tingkat yang baru secara kualitatif.



Dokumen Wolfram Mathematica asli dengan semua perhitungannya tersedia di GitHub .