Merancang filter Butterworth di Matlab dan mendapatkan koefisien [ab] filter sebagai bilangan bulat untuk generator kode Verilog HDL online

Saya telah merancang filter Butterworth low-pass yang sangat sederhana menggunakan Matlab. Cuplikan kode berikut menunjukkan apa yang telah saya lakukan.

fs = 2.1e6;
flow = 44 * 1000;
fNorm =  flow / (fs / 2);
[b,a] = butter(10, fNorm, 'low');

Dalam [b, a] disimpan koefisien filter. Saya ingin mendapatkan [b, a] sebagai bilangan bulat sehingga saya bisa menggunakan generator kode HDL online untuk menghasilkan kode di Verilog.

Nilai Matlab [b, a] tampaknya terlalu kecil untuk digunakan dengan pembuat kode online (skrip Perl sisi server menolak untuk menghasilkan kode dengan koefisien), dan saya bertanya-tanya apakah mungkin untuk mendapatkan [b, a] dalam bentuk yang dapat digunakan sebagai input yang tepat.

Koefisien yang saya dapatkan di Matlab adalah:

Koefisien b yang saya dapatkan di Matlab adalah:

1.0167e-012
1.0167e-011
4.5752e-011
1.2201e-010
2.1351e-010
2.5621e-010
2.1351e-010
1.2201e-010
4.5752e-011
1.0167e-011
1.0167e-012

Menggunakan generator online, saya ingin mendesain filter dengan bitwidth 12-bit dan bentuk filter I atau II. Saya tidak tahu apa yang dimaksud dengan "bit pecahan" pada tautan di atas.

Menjalankan generator kode (http://www.spiral.net/hardware/filter.html) dengan koefisien [b, a] yang tercantum di atas, dengan bit fraksional yang ditetapkan pada 20 dan bitwidth 12, saya menerima run error berikut :

Integer A constants: 1048576 -9603383 39613104 -96912015 155720456 -171714386 131597231 -69209161 23904282 -4896220 451621
Integer B constants: 0 0 0 0 0 0 0 0 0 0 0

Error: constants wider than 26 bits are not allowed, offending constant = -69209161, effective bitwidth = 7 mantissa + 20 fractional = 27 total.

An error has occurred - please revise the input parameters.

Bagaimana saya bisa mengubah desain saya sehingga kesalahan ini tidak terjadi?

UPDATE: Menggunakan Matlab untuk menghasilkan filter Butterworth urutan ke-6, saya mendapatkan koefisien berikut:

Untuk sebuah:

untuk b:

0.0064e-005
0.0382e-005
0.0954e-005
0.1272e-005
0.0954e-005
0.0382e-005
0.0064e-005

Menjalankan generator kode online (http://www.spiral.net/hardware/filter.html), saya sekarang menerima kesalahan berikut (dengan bit fraksional 8 dan bitwidth 20):

./iirGen.pl -A 256  '-1405' '3221' '-3943' '2719' '-1001' '153' -B  '0' '0' '0' '0' '0' '0' '0' -moduleName acm_filter -fractionalBits 8 -bitWidth 20 -inData inData  -inReg   -outReg  -outData outData -clk clk -reset reset -reset_edge negedge -filterForm 1  -debug  -outFile ../outputs/filter_1330617505.v 2>&1 
At least 1 non-zero-valued constant is required.  Please check the inputs and try again.

Mungkin koefisien b terlalu kecil, atau mungkin pembuat kode (http://www.spiral.net/hardware/filter.html) menginginkan [b, a] dalam format lain?

MEMPERBARUI:

Mungkin yang perlu saya lakukan adalah skala koefisien [b, a] dengan jumlah bit fraksional untuk mendapatkan koefisien sebagai bilangan bulat.

a .* 2^12
b .* 2^12

Namun, saya masih berpikir bahwa koefisien b sangat kecil. Apa yang saya lakukan salah di sini?

Mungkin tipe filter lain (atau metode desain filter) akan lebih cocok? Adakah yang bisa memberikan saran?

UPDATE: Seperti yang disarankan oleh Jason R dan Christopher Felton dalam komentar di bawah ini, filter SOS akan lebih cocok. Saya sekarang telah menulis beberapa kode Matlab untuk mendapatkan filter SOS.

fs = 2.1e6;
flow = 44 * 1000;      
fNorm =  flow / (fs / 2);
[A,B,C,D] = butter(10, fNorm, 'low');
[sos,g] = ss2sos(A,B,C,D);

Matriks SOS yang saya dapatkan adalah:

1.0000    3.4724    3.1253    1.0000   -1.7551    0.7705
1.0000    2.5057    1.9919    1.0000   -1.7751    0.7906
1.0000    1.6873    1.0267    1.0000   -1.8143    0.8301
1.0000    1.2550    0.5137    1.0000   -1.8712    0.8875
1.0000    1.0795    0.3046    1.0000   -1.9428    0.9598

Apakah mungkin untuk tetap menggunakan alat penghasil kode Verilog (http://www.spiral.net/hardware/filter.html) untuk mengimplementasikan filter SOS ini, atau haruskah saya menulis Verilog dengan tangan? Apakah tersedia referensi yang bagus?

Saya akan bertanya-tanya apakah filter FIR akan lebih baik untuk digunakan dalam situasi ini.

MOREOVER: Filter IIR rekursif dapat diimplementasikan menggunakan matematika integer dengan menyatakan koefisien sebagai fraksi. (Lihat buku pemrosesan sinyal DSP Smith yang sangat baik untuk perincian lebih lanjut: http://www.dspguide.com/ch19/5.htm )

Program Matlab berikut mengubah koefisien filter Butterworth menjadi bagian fraksional menggunakan fungsi Matlab rat (). Kemudian seperti yang disebutkan dalam komentar, bagian urutan kedua dapat digunakan untuk mengimplementasikan filter secara numerik (http://en.wikipedia.org/wiki/Digital_biquad_filter).

% variables
% variables
fs = 2.1e6;                     % sampling frequency           
flow = 44 * 1000;               % lowpass filter


% pre-calculations
fNorm =  flow / (fs / 2);       % normalized freq for lowpass filter

% uncomment this to look at the coefficients in fvtool
% compute [b,a] coefficients
% [b,a] = butter(7, fNorm, 'low');
% fvtool(b,a)  

% compute SOS coefficients (7th order filter)
[z,p,k] = butter(7, fNorm, 'low');

% NOTE that we might have to scale things to make sure
% that everything works out well (see zp2sos help for 'up' and 'inf' options)
sos = zp2sos(z,p,k, 'up', 'inf'); 
[n,d] = rat(sos); 
sos_check = n ./ d;  % this should be the same as SOS matrix

% by here, n is the numerator and d is the denominator coefficients
% as an example, write the the coefficients into a C code header file
% for prototyping the implementation

 % write the numerator and denominator matices into a file
[rownum, colnum] = size(n);  % d should be the same
sections = rownum;           % the number of sections is the same as the number of rows
fid = fopen('IIR_coeff.h', 'w');

fprintf(fid, '#ifndef IIR_COEFF_H\n');
fprintf(fid, '#define IIR_COEFF_H\n\n\n');
for i = 1:rownum
   for j = 1:colnum

       if(j <= 3)  % b coefficients
            bn = ['b' num2str(j-1) num2str(i) 'n' ' = ' num2str(n(i,j))];
            bd = ['b' num2str(j-1) num2str(i) 'd' ' = ' num2str(d(i,j))];
            fprintf(fid, 'const int32_t %s;\n', bn);
            fprintf(fid, 'const int32_t %s;\n', bd);

       end
       if(j >= 5)  % a coefficients
            if(j == 5) 
                colstr = '1'; 
            end
            if(j == 6) 
                colstr = '2'; 
            end
            an = ['a' colstr num2str(i) 'n' ' = ' num2str(n(i,j))];
            ad = ['a' colstr num2str(i) 'd' ' = ' num2str(d(i,j))];
            fprintf(fid, 'const int32_t %s;\n', an);
            fprintf(fid, 'const int32_t %s;\n', ad);
       end
   end
end

% write the end of the file
fprintf(fid, '\n\n\n#endif');
fclose(fid);

matlab Nicholas Kinar
sumber

Filter IIR tingkat tinggi seperti ini biasanya diimplementasikan menggunakan bagian urutan kedua ; Anda mendapatkan filter yang Anda inginkan dengan membuat beberapa tahapan orde kedua (dengan satu tahap orde satu jika urutan yang diinginkan aneh). Ini biasanya merupakan implementasi yang lebih kuat daripada langsung menerapkan filter tingkat tinggi.

Jason R

Jika Anda tidak melakukan apa yang disarankan @JasonR, Anda akan memiliki ukuran kata yang sangat besar. Filter seperti ini dapat gagal floating point presisi tunggal ketika diimplementasikan dengan struktur IIR dasar, Anda memerlukan SOS.

Christopher Felton

@JasonR: Terima kasih telah menyarankan ini. Saya telah memperbarui dengan jawaban di atas.

Nicholas Kinar

@ChristopherFelton: Terima kasih telah membantu memperkuat ini.

Nicholas Kinar

Ya dengan matriks SOS baru Anda, Anda dapat membuat 3 filter dari situs. Atau Anda dapat menggunakan kode saya di sini . Ini akan bekerja sama dengan situs web. Saya dengan senang hati akan memperbarui skrip ke kecuali matriks SOS.

Christopher Felton

Jawaban:

Seperti yang didiskusikan, yang terbaik adalah menggunakan penjumlahan bagian, yaitu memecah filter tingkat tinggi menjadi filter urutan kedua. Pertanyaan yang diperbarui memiliki matriks SOS. Menggunakan kode ini dan contoh di sini objek Python dapat digunakan untuk menghasilkan bagian individu.

Di matlab

save SOS

Dengan Python

import shutil
import numpy
from scipy.io import loadmat
from siir import SIIR

matfile = loadmat('SOS.mat')  
SOS = matfile['SOS']
b = numpy.zeros((3,3))
a = numpy.zeros((3,3))
section = [None for ii in range(3)]
for ii in xrange(3):
    b[ii] = SOS[ii,0:3]
    a[ii] = SOS[ii,3:6]

    section[ii] = SIIR(b=b[ii], a=a[ii], W=(24,0))
    section[ii].Convert()  # Create the Verilog for the section
    shutil.copyfile('siir_hdl.v', 'iir_sos_section%d.v'%(ii))

Informasi tambahan tentang titik tetap dapat ditemukan di sini

Christopher Felton
sumber

Terima kasih banyak untuk semua tautan penuh wawasan, dan untuk kode Python; Saya harap jawaban Anda (dan jawaban lain yang diposting di sini) berfungsi sebagai referensi yang baik untuk banyak orang lain. Saya berharap bahwa saya dapat menandai semua tanggapan di sini sebagai yang diterima.

Nicholas Kinar

Jika Anda memiliki masalah, beri tahu saya dan saya akan memperbarui / memperbaiki kode jika tidak berhasil untuk Anda. Saya akan memodifikasinya (relatif segera, doh) untuk langsung menerima matriks SOS.

Christopher Felton

Saya mencoba menerapkan versi saya sendiri dari contoh Anda. Di sistem saya, saya harus menggunakan "dari nol impor numpy" dan mengubah loatmat ke loadmat (). Apakah matriks SOS diberikan oleh Matlab ( mathworks.com/help/toolbox/signal/ref/ss2sos.html ) dalam format yang sama seperti yang diharapkan? Saya menerima kesalahan berikut ketika mencoba mengakses matriks SOS: "TypeError: tipe yang tidak dapat dihancurkan" ketika penerjemah mencapai garis "b [ii] = SOS [0: 3, ii]"

Nicholas Kinar

Itu akan tergantung pada format file SOS.mat. Jika Anda hanya >>> mencetak (matfile) itu akan menunjukkan kepada Anda kunci dalam file .mat yang dimuat. Scipy.io.loadmat selalu dimuat sebagai kamus (BOM).

Christopher Felton

Ya, itu benar, output dari 0 adalah input ke 1 dan seterusnya. Sedikit pemikiran perlu dimasukkan ke dalam kata lebar. Standarnya adalah dua fraksi menggunakan 24 bit (0 integer, 23 fraksi, 1 tanda). Saya yakin Anda awalnya ingin menggunakan lebar kata yang lebih kecil.

Christopher Felton

'Bit fraksional' adalah jumlah bit dalam bus yang telah Anda dedikasikan untuk mewakili bagian fraksional dari angka (misalnya, .75 dalam 3.75).

Katakanlah Anda memiliki lebar bus digital 4 bit, angka apa yang 1001diwakilinya? Ini bisa berarti '9' jika Anda memperlakukannya sebagai bilangan bulat positif (2 ^ 3 + 2 ^ 0 = 8 + 1 = 9). Atau bisa berarti -7 dalam notasi komplemen dua: (-2 ^ 3 + 2 ^ 0 = -8 + 1 = -7).

Bagaimana dengan angka dengan beberapa fraksi di dalamnya, yaitu angka 'nyata'? Bilangan real dapat direpresentasikan dalam perangkat keras sebagai "titik tetap" atau "titik mengambang". Sepertinya generator filter tersebut menggunakan titik tetap.

Kembali ke bus 4 bit kami (1001 ). Mari kita perkenalkan titik biner jadi kita dapatkan 1.001. Apa artinya ini adalah yang sekarang menggunakan bit pada RHS dari titik untuk membangun bilangan bulat, dan bit pada LHS untuk membangun sebagian kecil. Angka yang diwakili oleh bus digital yang disetel ke 1.0011.125 ( 1* 2 ^ 0 + 0* 2 ^ -1 + 0* 2 ^ -2 + 1* 2 ^ -3 = 1 + 0.125 = 1.125). Dalam hal ini, dari 4 bit dalam bus, kami menggunakan 3 bit untuk mewakili bagian fraksional dari angka. Atau, kami memiliki 3 bit fraksional.

Jadi, jika Anda memiliki daftar bilangan real seperti yang Anda miliki di atas, Anda sekarang harus memutuskan berapa banyak bit fraksional yang ingin Anda wakili. Dan inilah komprominya: semakin banyak bit fraksional yang Anda gunakan, semakin dekat Anda dapat mewakili angka yang Anda inginkan, tetapi semakin besar sirkuit Anda harus. Dan terlebih lagi, semakin sedikit bit fraksional yang Anda gunakan, semakin jauh respons frekuensi aktual filter akan menyimpang dari yang Anda rancang di awal!

Dan untuk memperburuk keadaan, Anda mencari untuk membangun filter Infinite Impulse Response (IIR). Ini sebenarnya bisa menjadi tidak stabil jika Anda tidak memiliki bit fraksional dan integer yang cukup!

Marty
sumber

Terima kasih telah memberikan jawaban mendalam ini. Saya sudah mencoba menjalankan pembuat kode menggunakan koefisien b kecil di atas, dan saya masih menerima beberapa kesalahan. Bisakah Anda menyarankan sesuatu yang bisa saya lakukan untuk menjalankan generator dengan benar? Saya akan memperbarui jawaban di atas untuk menunjukkan apa yang telah saya lakukan.

Jadi Marty telah mengurus dengan baik pertanyaan bit. Ke filter itu sendiri, saya pikir Anda kemungkinan mendapatkan peringatan atau keluhan dari matlab tentang koefisien yang berskala buruk? Ketika saya memplot filter, dari scipy bukan matlab tapi sepertinya sangat mirip.

Tanggapan

Yang 100 dB turun di passband! Jadi, Anda mungkin ingin memastikan Anda menginginkan filter pesanan yang lebih kecil, yang akan membantu implementasi Anda. Ketika saya sampai ke filter urutan ke-6 saya berhenti mendapatkan keluhan tentang koefisien buruk. Mungkin coba kurangi pesanan dan lihat apakah masih memenuhi persyaratan Anda.

Macduff
sumber

Terima kasih telah menyarankan ini! Saya pikir filter urutan ke-6 akan bekerja dengan baik. Menggunakan fvtool matlab, saya pikir responsnya bagus untuk aplikasi saya. Sekarang saya telah memperbarui respons saya di atas. Namun, ada yang salah dengan generator kode Verilog HDL ( spiral.net/hardware/filter.html ). Mungkin ia menginginkan [b, a] dalam format lain. Selain itu, +1 untuk penggunaan SciPy.