Bagaimana cara kerja filter low-pass secara programatik?

Saya telah mengerjakan filter low pass sederhana untuk pengukuran <100 Hz dalam aplikasi saya. Namun sejauh ini, saya berjuang dengan teori di balik itu semua. Sangat keren saya membuatnya bekerja, tetapi saya akan benar-benar menikmatinya jika saya tahu bagaimana / mengapa itu bekerja.

Saya menemukan kode berikut:

void getLPCoefficientsButterworth2Pole(const int samplerate, const double cutoff, double* const ax, double* const by)
{
    double PI = M_PI;
    double sqrt2 = sqrt(2);

    double QcRaw  = (2 * PI * cutoff) / samplerate; // Find cutoff frequency in [0..PI]
    double QcWarp = tan(QcRaw); // Warp cutoff frequency

    double gain = 1 / ( 1 + sqrt2 / QcWarp + 2 / ( QcWarp * QcWarp ) );

    by[2] = ( 1 - sqrt2 / QcWarp + 2 / ( QcWarp * QcWarp ) ) * gain;
    by[1] = ( 2 - 2 * 2 / ( QcWarp * QcWarp ) ) * gain;
    by[0] = 1;

    ax[0] = 1 * gain;
    ax[1] = 2 * gain;
    ax[2] = 1 * gain;
}

Untuk menghitung Koefisien. Kemudian, dalam sampel audio, saya 'lulus rendah' ​​dengan cara ini:

        xv[2] = xv[1];
        xv[1] = xv[0];

        xv[0] = pData[j];
        yv[2] = yv[1];
        yv[1] = yv[0];

        yv[0] = (ax[0] * xv[0] + ax[1] * xv[1] + ax[2] * xv[2]
                   - by[1] * yv[0]
                   - by[2] * yv[1]);

        pData[j] = yv[0];

Untuk mendapatkan desain lowpass.

Saya bertanya-tanya beberapa hal:

1. Saya menerima sampel audio dalam float * array sederhana. Berapa nomor float itu? Satu-satunya hal yang saya lihat adalah angka, bagaimana itu sepotong suara?
2. Kode menggunakan perhitungan sebelumnya (tiga di antaranya) dalam perhitungan baru per sampel. Apakah itu berarti bahwa 2 sampel data pertama tidak difilter dengan benar? (Bukan berarti itu penting karena hanya 2 sampel, tetapi hanya ingin tahu)
3. Mencoba mempelajari semuanya, saya menemukan beberapa formula untuk filter Butterworth (2nd Pole). Bagaimana rumus-rumus itu tercermin dalam kode ini? Tidak ada rumus yang saya temukan memiliki perhitungan ini yang dapat Anda lihat di fungsi 'getLPCoefficientsButterworth2Pole ()'.

1. Nomor float * array adalah penunjuk ke array. Ini adalah nomor tunggal yang berisi alamat elemen pertama dari array nilai float.

2. Biasanya, kondisi awal (yaitu elemen 'masa lalu' awal dari x dan y) adalah 0, tetapi jika nilainya tidak sama dengan 0 maka tidak ada masalah besar juga, karena setelah beberapa saat kondisi awal tidak berpengaruh pada output sinyal untuk setiap filter stabil. Dan filter Anda jelas stabil.

3. Fungsi transfer lowpass urutan kedua dengan karakteristik Butterworth dan frekuensi cut-off (dengan dalam Hertz) diberikan oleh$\omega_a=2*\pi*f_a$$f_a$

$$H(s) = \frac{\omega_a^2}{s^2+\sqrt{2}\omega_a s + \omega_a^2}\tag{1}$$ Ini adalah hasil standar yang dapat Anda temukan dengan mudah di web . Untuk mendapatkan filter waktu diskrit, Anda dapat menerapkan apa yang disebut transformasi bilinear:

$$s = 2f_s\frac{z-1}{z+1}\tag{2}$$

di mana adalah frekuensi sampling. Ini diperlukan karena sumbu frekuensi dari sinyal analog ( ) perlu dipetakan ke rentang frekuensi yang diizinkan untuk sinyal waktu diskrit ( ). Karena transformasi ini mengurangi frekuensi, kita perlu menghitung frekuensi cut-off analog yang diinginkan dari frekuensi cut-off yang diberikan dalam domain waktu diskrit:$f_s$$0\le f\le\infty$$0\le f\le f_s/2$$f_d$

$$\omega_a = 2f_s \tan(\frac{\omega_d}{2})\text{ with }\omega_d = 2\pi f_d/f_s$$

Jika Anda memasukkan (2) ke dalam (1) Anda dapatkan

$$H(z) = k\frac{z^2+b_1z+b_2}{z^2+a_1z+a_2}\tag{3}$$

dengan

$$k=\frac{\alpha^2}{1+\sqrt{2}\alpha+\alpha^2}\text{ (gain)}$$ $$a_1 = \frac{2(\alpha^2-1)}{1+\sqrt{2}\alpha+\alpha^2},\quad a_2=\frac{1-\sqrt{2}\alpha+\alpha^2}{1+\sqrt{2}\alpha+\alpha^2},\quad b_1=2,\quad b_2=1$$

tempat saya menggunakan . Fungsi transfer (3) sesuai dengan persamaan domain waktu$\alpha=\tan(\frac{\omega_d}{2})$

$$y(n) = kx(n) + kb_1x(n-1) + kb_2x(n-2) - a_1y(n-1) - a_2y(n-2)$$

Seperti yang Anda lihat, ada beberapa kemiripan dengan kode Anda. Namun, ada beberapa perbedaan juga dan Anda harus memeriksa respons frekuensi filter Anda. Saya pikir persamaan di atas benar, tetapi terserah Anda untuk memeriksanya dan memutuskan versi mana yang benar.

Anda bertanya bagaimana filter low pass bekerja dan menyebutkan bahwa filter menggunakan nilai sebelumnya dari data Anda. Ini adalah diskusi non-teknis tentang apa yang terjadi pada filter low pass.

Filter low pass mengambil tampilan berbeda (bergeser dalam waktu) dari sinyal Anda, menskala mereka dan menambahkannya bersama-sama. Anda dapat membayangkan menggambar sinyal Anda 3 kali, satu menjadi arus, yang kedua digeser oleh satu waktu sampel, yang ketiga digeser oleh 2 kali sampel.

Pada frekuensi rendah, semua tampilan tampak sangat mirip (bergeser oleh sampel tunggal nyaris tidak mengubah di mana Anda berada pada sinyal kapan saja dalam waktu singkat). Dalam hal ini, ketiga versi akan ditambahkan bersama dengan cara yang konstruktif (atau setidaknya tidak merusak), sehingga sinyal melewati filter.

Sekarang pindah ke frekuensi yang lebih tinggi, setiap versi sinyal yang bergeser menjadi lebih berbeda pada setiap saat (titik sampel) tertentu dan bahkan mungkin terbalik dalam tanda. Pada frekuensi yang lebih tinggi ini, tiga versi sinyal Anda cenderung untuk membatalkan (ditambahkan secara destruktif) sehingga sinyal menjadi dilemahkan.

Berbagai jenis filter mengatur jenis gangguan konstruktif / destruktif ini terjadi pada pita frekuensi yang sesuai untuk membuat filter low-pass, band-pass atau high-pass.

Itu tergantung pada bagaimana Anda menginginkannya. Bagi saya, saya telah menerapkan di C. Jadi, apa yang saya lakukan adalah saya membuat koefisien filter dalam matlabmenggunakan firls. Lalu saya menyimpan koefisien filter dalam array dan kemudian meneruskan koefisien filter ini ke fungsi filter dalam C. Di matlabdalamnya cukup mudah untuk menghasilkan koefisien filter, dan kemudian menggunakan operasi konvolusi atau menggunakan filter tertentu. Tetapi di dalam Ckamu harus menerapkan filter juga. Untuk bagian saya, saya harus melakukannya untuk DSP, jadi diimplementasikan dalam Cdan mendapatkan hasil dalam matlabmenggunakan fungsi mex. Untuk menghitung koefisien, gunakan perintah ini di matlab:

n=17                  %filter order
f=[0 0.167 0.333 1]   %Frequency band edges
a=[1 1 0 0]           %Desired amplitudes
fir= firls (n, f,a )