Cara yang Benar untuk Melakukan Sinc Downsampling (DFT Downsampling) untuk Sinyal Diskrit Sampel Berseragam dengan Jumlah Sampel Terbatas
Diberikan sinyal $ \left\{ x [ 0 ], x [ 1 ], ..., x [ N - 1 ] \right\} $ apa cara yang benar untuk downsample dalam domain frekuensi (interpolasi Sinc)?
Jawaban
Interpolasi dalam Frekuensi (Domain DFT)
Penerapannya sudah terkenal. Di MATLAB akan menjadi seperti ini:
if(numSamplesO > numSamples)
% Upsample
halfNSamples = numSamples / 2;
if(mod(numSamples, 2) ~= 0) % Odd number of samples
vXDftInt = interpFactor * [vXDft(1:ceil(halfNSamples)); zeros(numSamplesO - numSamples, 1, 'like', vXDft); vXDft((ceil(halfNSamples) + 1):numSamples)];
else % Even number of samples -> Special Case
vXDftInt = interpFactor * [vXDft(1:halfNSamples); vXDft(halfNSamples + 1) / 2; zeros(numSamplesO - numSamples - 1, 1, 'like', vXDft); vXDft(halfNSamples + 1) / 2; vXDft((halfNSamples + 2):numSamples)];
end
else
% Downsample
halfNSamples = numSamplesO / 2;
if(mod(numSamples, 2) ~= 0) % Odd number of samples
vXDftInt = interpFactor * [vXDft(1:ceil(halfNSamples)); vXDft((numSamples - floor(halfNSamples) + 1):numSamples)];
else % Even number of samples -> Special Case
vXDftInt = interpFactor * [vXDft(1:halfNSamples); vXDft(halfNSamples + 1) / 2; vXDft((numSamples - halfNSamples + 2):numSamples)];
end
end
Jadi kami menangani 2 kasus di sini:
- Upsample
Kami menambahkan sampel nol ke bagian tengah DFT untuk mencocokkan jumlah sampel output (numSamplesO).
Kami menangani kasus jumlah input sampel (numSamples) genap. Dalam hal ini kami membagi sampel Nyquist ($ X \left[ N / 2 \right] $) menjadi 2 di mana $ N $ adalah jumlah sampel yang diinput. - Downsample
Kami menghapus sampel dari bagian tengah DFT untuk mencocokkan jumlah sampel output (numSamplesO).
Kami menangani kasus jumlah keluaran sampel (numSamplesO) genap. Dalam hal ini kami membagi sampel menjadi Nyquist ($ X \left[ M / 2 \right] $) menjadi 2 di mana $ M $ adalah jumlah keluaran sampel.
Pertanyaannya adalah, mengapa kita melakukannya seperti ini? Mengapa faktor interpolasi interpFactor? Dimana faktor pemisahan$ 0.5 $berasal dari?
Untuk menjawab itu perlu kita ingat DFT pada dasarnya adalah Discrete Fourier Series (DFS).
Artinya asumsi yang paling penting adalah data bersifat periodik baik dalam domain waktu maupun frekuensi.
Sekarang, karena DFT pada dasarnya adalah DFS , cara alami untuk menginterpolasi sinyal dalam periode tersebut akan menggunakan Seri Fourier.
Sebelum membahas detailnya, mari kita tentukan 2 set bilangan bulat yang akan digunakan untuk menentukan nilai indeks:
$$ \begin{aligned} \mathcal{K}_{DFS}^{N} & = \left\{- \left\lceil \frac{N - 1}{2} \right\rceil, - \left\lceil \frac{N - 1}{2} \right\rceil + 1, \ldots, -1, 0, 1, \ldots, \left\lceil \frac{N - 1}{2} \right\rceil - 1, \left\lceil \frac{N - 1}{2} \right\rceil \right\} \\ \mathcal{K}_{DFT}^{N} & = \left\{- \left\lceil \frac{N - 1}{2} \right\rceil, - \left\lceil \frac{N - 1}{2} \right\rceil + 1, \ldots, -1, 0, 1, \ldots, \left\lceil \frac{N - 1}{2} \right\rceil - 1, \left\lfloor \frac{N - 1}{2} \right\rfloor \right\} \\ \end{aligned} $$
Artinya, untuk sinyal dengan bandwidth maksimum $ \frac{1}{2 T} $ diambil sampelnya dengan Teorema Sampling untuk $ t \in \left[ 0, N T \right) $ dimana $ T $ adalah periode pengambilan sampel dan $ P = N T $ adalah periode fungsi:
$$ \begin{aligned} x \left( t \right) {\Big|}_{t = n T} & = \sum_{k = - \left\lceil \frac{N - 1}{2} \right\rceil}^{\left\lceil \frac{N - 1}{2} \right\rceil} {c}_{k} {e}^{ j 2 \pi \frac{k t}{P} } && \text{By Fourier Series} \\ & = \sum_{k = - \left\lceil \frac{N - 1}{2} \right\rceil}^{\left\lceil \frac{N - 1}{2} \right\rceil} {c}_{k} {e}^{ j 2 \pi \frac{k t}{N T} } && \text{By the period of the function / series} \\ & = \sum_{k = - \left\lceil \frac{N - 1}{2} \right\rceil}^{\left\lceil \frac{N - 1}{2} \right\rceil} {c}_{k} {e}^{ j 2 \pi \frac{k n}{N} } && \text{Setting $ t = n T $} \\ & = \frac{1}{N} \sum_{k = - \left\lceil \frac{N - 1}{2} \right\rceil}^{\left\lfloor \frac{N - 1}{2} \right\rfloor} X \left[ k \right] {e}^{ j 2 \pi \frac{k n}{N} } && \text{The DFT} \end{aligned} $$
Rumus di atas berfungsi untuk kasus genap $ N = 2 l, \; l \in \mathbb{N} $ dan untuk kasus aneh $ N = 2 l + 1, \; l \in \mathbb{N} $. Di atas mendefinisikan hubungan antara koefisien DFT dan Koefisien Deret Fourier :
$$ {c}_{k} = \begin{cases} \frac{ X \left[ k \right ] }{2 N} & \text{ if } k = \frac{N}{2} \\ \frac{ X \left[ k \right ] }{2 N} & \text{ if } k = -\frac{N}{2} \\ \frac{ X \left[ k \right ] }{N} & \text{ if } k \notin \left\{\frac{N}{2}, -\frac{N}{2} \right\} \end{cases}, \; k \in \mathcal{K}_{DFS}^{N} $$
Tetapi tidak ada yang menghentikan kami untuk menggunakan titik pengambilan sampel lain untuk set mana pun $ { \left\{ {t}_{m} \right\}}_{m = 0}^{M - 1} $ dimana $ \forall m, {t}_{m} \in \left[ 0, N T \right) $. Yang memberikan$ x \left( t \right) = \frac{1}{N} \sum_{k = - \left\lceil \frac{N - 1}{2} \right\rceil}^{\left\lfloor \frac{N - 1}{2} \right\rfloor} X \left[ k \right] {e}^{ j 2 \pi \frac{k t}{N T} } $ untuk $ t \in \left[ 0, N T \right) $. Ini akan berfungsi untuk sinyal yang kompleks dan nyata.
Untuk sinyal nyata,$ x \left( t \right) \in \mathbb{R} $kita juga bisa menggunakan bentuk Cosine dari DFT :
$$ \begin{aligned} x \left( t \right) & = \sum_{k = - \left\lceil \frac{N - 1}{2} \right\rceil}^{\left\lceil \frac{N - 1}{2} \right\rceil} {c}_{k} {e}^{ j 2 \pi \frac{k t}{N T} } && \text{From the above} \\ & = \sum_{k = - \left\lceil \frac{N - 1}{2} \right\rceil}^{\left\lceil \frac{N - 1}{2} \right\rceil} \left| {c}_{k} \right| \cos \left( 2 \pi \frac{k t}{N T} + \angle {c}_{k} \right) && \text{Fourier series in its Cosine form} \\ & = \sum_{k = - \left\lceil \frac{N - 1}{2} \right\rceil}^{\left\lfloor \frac{N - 1}{2} \right\rfloor} \frac{\left| X \left[ k \right] \right|}{N} \cos \left( 2 \pi \frac{k t}{N T} + \angle X \left[ k \right] \right) && \text{Fourier series in its Cosine form} \\ & = \sum_{k = 0}^{\left\lfloor \frac{N - 1}{2} \right\rfloor} {\alpha}_{k} \frac{\left| X \left[ k \right] \right|}{N} \cos \left( 2 \pi \frac{k t}{N T} + \angle X \left[ k \right] \right) && \text{Using the DFT conjugate symmetry of a real signal} \end{aligned} $$
Dimana $ {\alpha}_{k} = \begin{cases} 1 & \text{ if } k \in \left\{ 0, \frac{N}{2} \right\} \\ 2 & \text{ else } \end{cases} $.
Jadi sekarang kita perlu memikirkan apa yang kita lihat di sini dan bagaimana hubungannya dengan algoritma di atas.
Pertama-tama kita perlu memperhatikan bahwa trik utama di sini adalah bahwa bentuk asli DFT harus saat indeks berjalan$ k \in \mathcal{K}_{DFT}^{N} $. Maka lebih mudah untuk melihat koneksi ke asal-usul Discrete Fourier Series ( DFS ) dari DFT .
Catatan : Dalam praktiknya, DFT didefinisikan (Dan dihitung) dengan$ k \in \left\{ 0, 1, \ldots, N - 1 \right\} $.
Jika kita memilih himpunan grid waktu seragam keluaran $ { \left\{ {t}_{m} \right\}}_{m = 0}^{M - 1} $ untuk menjadi dalam bentuk $ {t}_{m} = m {T}_{s} $ di mana tingkat upsampling (Kami akan menangani downsampling nanti) $ q = \frac{M}{N} \geq 1 $maka jelas apa yang perlu dilakukan dengan melihat IDFT untuk memulihkan grid:
$$ x \left[ m \right] = \frac{1}{M} \sum_{k = 0}^{M - 1} \tilde{X} \left[ k \right] {e}^{j 2 \pi \frac{k m}{M}} = \frac{1}{M} \sum_{k = - \left\lceil \frac{M - 1}{2} \right\rceil}^{\left\lfloor \frac{M - 1}{2} \right\rfloor} \tilde{X} \left[ k \right] {e}^{j 2 \pi \frac{k m}{M}} $$
Sekarang kita perlu mencocokkan ini dengan rumus interpolasi dari atas. Karena itu adalah transformasi linier dikalikan dengan$ q $akan menjaga konstanta. Kami juga bisa memperhatikan itu$ \forall m, \frac{m}{M} = \frac{{t}_{m}}{N T} $ karenanya dengan mengatur:
$$ \tilde{X} \left[ k \right] = \begin{cases} X \left[ k \right] & \text{ if } k \in \mathcal{K}_{DFT}^{N} \setminus \left\{ k \mid k = \frac{N}{2} \right\} \\ \frac{X \left[ k \right]}{2} & \text{ if } k = \frac{N}{2} \\ 0 & \text{ if } k \notin \mathcal{K}_{DFT}^{N} \end{cases} $$
Dari $ N $ periodisitas dari DFT kita dapat menulis interpolasi akhir untuk grid waktu yang seragam dengan faktor interpolasi $ q $:
$$ x \left[ m \right] = \frac{q}{M} \sum_{k = 0}^{M - 1} \hat{X} \left[ k \right] {e}^{j 2 \pi \frac{k m}{M}} $$
Dimana $ \hat{X} \left[ k \right] $ didefinisikan sebagai:
$$ \hat{X} \left[ k \right] = \begin{cases} X \left[ k \right] & \text{ if } k \in \left\{ 0, 1, \ldots, N - 1 \right\} \setminus \left\{ \frac{N}{2} \right\} \\ \frac{X \left[ k \right]}{2} & \text{ if } k = \frac{N}{2} \\ 0 & \text{ if } k \in \left\{ N, N + 1, \ldots, M - 1 \right\} \end{cases} $$
Yang persis seperti yang kami lakukan pada kode upsample di atas.
Bagaimana dengan downsample? Nah, kita bisa menggunakan intuisi yang sama di domain DFT seperti yang ditunjukkan oleh kode. Hal ini pada dasarnya karena interpolasi menggunakan koefisien Deret Fourier tidak lain adalah perkalian dalam domain frekuensi dengan Kernel Dirichlet yang merupakan ekuivalen periodik dari$ \operatorname{sinc} \left( \cdot \right) $fungsi. Ini juga merupakan intuisi untuk$ \frac{1}{2} $faktor, karena kita mengalikan dengan persegi dengan nilai 1 pada domain frekuensi yang memiliki diskontinuitas lompat . Seri Fourier memang menyatu dengan nilai rata-rata lompatan saat penghentian. Sejak kita pergi dari$ 1 $ untuk $ 0 $, artinya nilai pada lompatan tersebut $ 0.5 $.
Jadi kode downsmaplign dan upsampling di atas hanya menerapkan Kernel Dirichlet ke data sesuai dengan frekuensi sampling dari input, pada kasus upsample dan output pada kasus downsample.
Metode lain untuk downsample adalah upampling ke faktor integer dari jumlah keluaran sampel. Kemudian gunakan desimasi (Ambil setiap ... sampel) untuk mendapatkan sampel. Angka 2 akan cocok untuk kasus data tidak memiliki energi dalam frekuensi antara laju rendah dan laju sampel. Jika ya, mereka tidak akan cocok.
Saya akan menambahkan Kode MATLAB ...
Catatan : Jawaban ini juga mencakup Upsampling . Harap pertimbangkan untuk membuka pertanyaan lain tentang Upsampling atau perluas pertanyaan ini.