Entwerfen des Butterworth-Filters in Matlab und Abrufen der Filter [ab] -Koeffizienten als Ganzzahlen für den Online-Verilog-HDL-Codegenerator

15

Ich habe mit Matlab einen sehr einfachen Butterworth-Tiefpassfilter entwickelt. Der folgende Codeausschnitt zeigt, was ich getan habe.

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

In [b, a] sind die Filterkoeffizienten gespeichert. Ich möchte [b, a] als Ganzzahl erhalten, damit ich mit einem Online-HDL-Codegenerator Code in Verilog generieren kann.

Die Matlab-Werte [b, a] scheinen zu klein für den Online-Codegenerator zu sein (das serverseitige Perl-Skript weigert sich, Code mit den Koeffizienten zu generieren), und ich frage mich, ob es möglich wäre, [b] zu erhalten. a] in einer Form, die als richtige Eingabe verwendet werden kann.

Die a-Koeffizienten, die ich in Matlab erhalte, sind:

1.0000
-9.1585
37.7780
-92.4225
148.5066
-163.7596
125.5009
-66.0030
22.7969
-4.6694
0.4307

Die b-Koeffizienten, die ich in Matlab erhalte, sind:

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

Mit dem Online-Generator möchte ich einen Filter mit einer 12-Bit-Bit-Breite und einer I- oder II-Filterform entwerfen. Ich weiß nicht, was unter "Nachkommastellen" im obigen Link zu verstehen ist.

Wenn ich den Codegenerator (http://www.spiral.net/hardware/filter.html) mit den oben aufgelisteten Koeffizienten [b, a] mit Bruchbits von 20 und einer Bitbreite von 12 ausführe, erhalte ich den folgenden Ausführungsfehler :

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. 

Wie kann ich mein Design ändern, damit dieser Fehler nicht auftritt?

UPDATE: Mit Matlab zur Erzeugung eines Butterworth-Filters 6. Ordnung erhalte ich folgende Koeffizienten:

Für ein:

1.0000
-5.4914
12.5848
-15.4051
10.6225
-3.9118
0.6010 

für b:

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

Beim Ausführen des Online-Codegenerators (http://www.spiral.net/hardware/filter.html) wird jetzt die folgende Fehlermeldung angezeigt (mit Nachkommastellen von 8 und einer Bitbreite von 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.

Vielleicht sind die b-Koeffizienten zu klein, oder vielleicht will der Codegenerator (http://www.spiral.net/hardware/filter.html) das [b, a] in einem anderen Format?

AKTUALISIEREN:

Vielleicht muss ich die [b, a] -Koeffizienten durch die Anzahl der Bruchbits skalieren, um die Koeffizienten als ganze Zahlen zu erhalten.

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

Ich denke jedoch immer noch, dass die b-Koeffizienten extrem klein sind. Was mache ich hier falsch?

Vielleicht ist ein anderer Filtertyp (oder eine andere Filterdesignmethode) besser geeignet? Könnte jemand einen Vorschlag machen?

UPDATE: Wie von Jason R und Christopher Felton in den Kommentaren unten vorgeschlagen, wäre ein SOS-Filter besser geeignet. Ich habe jetzt Matlab-Code geschrieben, um einen SOS-Filter zu erhalten.

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);

Die SOS-Matrix, die ich bekomme, ist:

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

Ist es weiterhin möglich, das Verilog-Codegenerierungstool (http://www.spiral.net/hardware/filter.html) zum Implementieren dieses SOS-Filters zu verwenden, oder sollte ich den Verilog einfach von Hand schreiben? Ist eine gute Referenz vorhanden?

Ich würde mich fragen, ob ein FIR-Filter in dieser Situation besser zu verwenden wäre.

WEITERES: Rekursive IIR-Filter können mit ganzzahliger Mathematik implementiert werden, indem Koeffizienten als Brüche ausgedrückt werden. (Weitere Informationen finden Sie in Smiths ausgezeichnetem DSP-Signalverarbeitungsbuch: http://www.dspguide.com/ch19/5.htm )

Das folgende Matlab-Programm konvertiert Butterworth-Filterkoeffizienten mit der Funktion Matlab rat () in Bruchteile. Wie in den Kommentaren erwähnt, können Abschnitte zweiter Ordnung verwendet werden, um den Filter numerisch zu implementieren (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);
Nicholas Kinar
quelle
4
IIR-Filter höherer Ordnung wie dieses werden normalerweise unter Verwendung von Abschnitten zweiter Ordnung implementiert ; Sie erhalten den gewünschten Filter, indem Sie mehrere Stufen zweiter Ordnung hintereinander schalten (mit einer einzigen Stufe erster Ordnung, wenn die gewünschte Reihenfolge ungerade ist). Dies ist normalerweise eine robustere Implementierung als die direkte Implementierung des Filters höherer Ordnung.
Jason R
3
Wenn Sie nicht das tun, was @JasonR vorschlägt, werden Sie sehr große Wörter haben. Filter wie diese können Gleitkommazahlen mit einfacher Genauigkeit verfehlen, wenn sie mit einer grundlegenden IIR-Struktur implementiert werden. Sie benötigen das SOS.
Christopher Felton
@JasonR: Danke, dass du dies vorschlägst. Ich habe per Antwort oben aktualisiert.
Nicholas Kinar
@ChristopherFelton: Vielen Dank, dass Sie dazu beigetragen haben, dies zu verstärken.
Nicholas Kinar
Ja, mit Ihrer neuen SOS-Matrix können Sie 3 Filter auf der Website erstellen. Oder Sie können meinen Code hier verwenden . Es funktioniert genauso wie die Website. Gerne aktualisiere ich das Skript bis auf SOS-Matrix.
Christopher Felton

Antworten:

5

Wie bereits erwähnt, ist es am besten, die Summe der Abschnitte zu verwenden, dh den Filter höherer Ordnung in kaskadierte Filter zweiter Ordnung zu unterteilen. Die aktualisierte Frage hat die SOS-Matrix. Mit diesem Code und einem Beispiel hier kann das Python-Objekt verwendet werden, um die einzelnen Abschnitte zu generieren.

In Matlab

save SOS

In 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))

Weitere Informationen zum Fixpunkt finden Sie hier

Christopher Felton
quelle
Vielen Dank für all die aufschlussreichen Links und für den Python-Code. Ich hoffe, dass Ihre Antwort (und die anderen hier veröffentlichten Antworten) vielen anderen als Referenz dienen. Ich wünschte, ich könnte alle Antworten hier als angenommen markieren.
Nicholas Kinar
1
Wenn Sie irgendwelche Probleme haben, lassen Sie es mich wissen und ich werde den Code aktualisieren / korrigieren, wenn es bei Ihnen nicht funktioniert. Ich werde es ändern (relativ bald, doh), um eine SOS-Matrix direkt zu akzeptieren.
Christopher Felton
1
Ich habe versucht, meine eigene Version aus Ihrem Beispiel zu implementieren. Auf meinem System musste ich "from numpy import zero" verwenden und loatmat in loadmat () ändern. Hat die von Matlab ( mathworks.com/help/toolbox/signal/ref/ss2sos.html ) angegebene SOS-Matrix dasselbe Format wie erwartet? Ich erhalte die folgende Fehlermeldung, wenn ich versuche, auf die SOS-Matrix zuzugreifen: "TypeError: Unhashable Type", wenn der Interpreter die Zeile "b [ii] = SOS [0: 3, ii]" erreicht
Nicholas Kinar
1
Das würde vom Format der SOS.mat-Datei abhängen. Wenn Sie einfach >>> drucken (matfile), werden Ihnen die Schlüssel in der geladenen .mat-Datei angezeigt. Die scipy.io.loadmat wird immer als Wörterbuch (BOMK) geladen.
Christopher Felton
1
Ja, das ist richtig, die Ausgabe von 0 ist die Eingabe von 1 und so weiter. Ein kleiner Gedanke muss in die Wortbreite gesteckt werden. Die Standardeinstellung ist 24-Bit-Bruch (0 Ganzzahl, 23 Bruch, 1 Vorzeichen). Ich glaube, Sie wollten ursprünglich eine kleinere Wortbreite verwenden.
Christopher Felton
10

Die 'Bruchbits' sind die Anzahl der Bits in einem Bus, die Sie für den Bruchteil einer Zahl reserviert haben (z. B. die 0,75 in 3,75).

Angenommen, Sie haben einen digitalen Bus, der 4 Bit breit ist. Welche Zahl steht dafür 1001? Es könnte '9' bedeuten, wenn Sie es als positive Ganzzahl behandeln (2 ^ 3 + 2 ^ 0 = 8 + 1 = 9). Oder es könnte -7 in der Zweierkomplementnotation bedeuten: (-2 ^ 3 + 2 ^ 0 = -8 + 1 = -7).

Was ist mit Zahlen mit einigen Brüchen, dh "reellen" Zahlen? Reelle Zahlen können in der Hardware als "Festkomma" oder "Gleitkomma" dargestellt werden. Es sieht so aus, als ob diese Filtergeneratoren Festpunkte verwenden.

Zurück zu unserem 4 Bit Bus ( 1001). Lassen Sie uns einen Binärpunkt einführen, damit wir erhalten 1.001. Dies bedeutet, dass jetzt die Bits auf der rechten Seite des Punkts verwendet wurden, um Ganzzahlen zu erstellen, und die Bits auf der linken Seite, um einen Bruch zu erstellen. Die Zahl, die durch einen Digitalbus dargestellt wird, der auf gesetzt ist, 1.001ist 1,125 ( 1* 2 ^ 0 + 0* 2 ^ -1 + 0* 2 ^ -2 + 1* 2 ^ -3 = 1 + 0,125 = 1,125). In diesem Fall verwenden wir von den 4 Bits im Bus 3, um den Bruchteil einer Zahl darzustellen. Oder wir haben 3 gebrochene Bits.

Wenn Sie also eine Liste mit reellen Zahlen wie oben haben, müssen Sie jetzt entscheiden, wie viele Bruchbits Sie darstellen möchten. Und hier ist der Kompromiss: Je mehr Bruchbits Sie verwenden, desto genauer können Sie die gewünschte Zahl darstellen, aber desto größer muss Ihre Schaltung sein. Und je weniger Bruchbits Sie verwenden, desto stärker weicht der tatsächliche Frequenzgang des Filters von dem ab, den Sie zu Beginn entworfen haben!

Und um die Sache noch schlimmer zu machen, möchten Sie einen Infinite Impulse Response (IIR) -Filter bauen. Diese können tatsächlich instabil werden, wenn Sie nicht genügend gebrochene und ganzzahlige Bits haben!

Marty
quelle
Vielen Dank für diese aufschlussreiche Antwort. Ich habe versucht, den Codegenerator mit den kleinen b-Koeffizienten oben auszuführen, und es werden immer noch einige Fehler angezeigt. Könnten Sie mir etwas vorschlagen, das ich tun könnte, um den Generator richtig laufen zu lassen? Ich werde die Antwort oben aktualisieren, um zu zeigen, was ich getan habe.
10

Marty hat sich also gut um die Frage der Bits gekümmert. Ich glaube, Sie erhalten auf dem Filter selbst eine Warnung oder Beschwerde von matlab über schlecht skalierte Koeffizienten? Wenn ich den Filter zeichne, von scipy nicht matlab, aber es ist wahrscheinlich sehr ähnlich.

Antwort

Das sind 100 dB im Durchlassbereich! Sie können also sicherstellen, dass Sie einen Filter mit kleinerer Reihenfolge wünschen, der Ihnen trotzdem bei der Implementierung hilft. Wenn ich zu einem Filter 6. Ordnung komme, höre ich auf, mich über schlechte Koeffizienten zu beschweren. Versuchen Sie möglicherweise, die Reihenfolge zu verringern, und prüfen Sie, ob sie noch Ihren Anforderungen entspricht.

macduff
quelle
Vielen Dank für den Vorschlag! Ich denke, dass ein Filter 6. Ordnung genauso gut funktionieren würde. Mit dem fvtool von matlab finde ich, dass die Antwort für meine Anwendung gut ist. Ich habe jetzt meine Antwort oben aktualisiert. Mit dem Verilog HDL-Codegenerator ( spiral.net/hardware/filter.html ) läuft jedoch immer noch etwas schief . Vielleicht will es das [b, a] in einem anderen Format. Zusätzlich +1 für die Verwendung von SciPy.