IIR Filter 2. Ordnung in FPGA implementieren

[Marcel D. (diablokiller999)]

Hi Leutz!
Nachdem ich erfolgreich ein FIR Filter in den FPGA gebracht habe, wollte 
ich nun schauen, ob ich mit einem IIR Filter Platz spare. Dazu habe ich 
mir ein IIR Filter mit Chebyshev 2 Eigenschaften in Matlab erstellt, 
Zähler/Nennerpolynom mithilfe der tf2sos-Funktion in ein SOS Filter 
überführt, der den Gain direkt zu den Koeffizienten multipliziert.

Zur Implementierung auf dem FPGA in einem getakteten Prozess habe ich 
mir nun überlegt, zuerst die Signale mit den Koeffizienten zu 
multiplizieren (nA1,nA2,nB1,nB2), diese dann zu summieren 
(nSUMA1A2,nSUMB1B2), die Summe mit dem Eingangssignal zu bilden (nSUMX), 
die beiden Speicherregister zu aktualisieren und das Produkt aus 
Eingangssumme und Koeffizient B0 zu berechnen sowie letztlich das 
Ausgangssignal aus der Summe von nB0 und nSUMB1B2 zu erhalten (siehe 
Bild). Da ich Samples mit höchstens 1MHz bekomme und mein FPGA mit 
100MHz läuft, habe ich also genügend Zeit alles in einem getakteten 
Prozess zu verarbeiten.

Nur leider sieht die Simulation kein bisschen nach einem IIR Filter aus, 
gebe ich einen Sprung ins System, erhalte ich nicht das für ein IIR 
Filter typische(?) Einschwingen. Wahrscheinlich ist schon mein Ansatz 
falsch gewählt? Momentan nutze ich zwei Koeffizienten, später soll das 
Design es eigentlich ermöglichen auch ein Filter 4. oder 6. Ordnung 
durch Kaskadierung zu erreichen (SOS Filter halt).

[Martin O. (ossi-2)]

Zeig mal deinen Code zusammen mit den Koeffizienten. Bei IIR Filtern ist 
es wichtig die "internen Werte" passend zu skalieren, damit alles 
klappt.

[Marcel D. (diablokiller999)]

Stehen in der angehängten IIR.vhd inklusive der Umrechnung in 32Bit 
signed Werte als Kommentar.

B0 : integer := 14419;         -- = ((2^31)/1,995) * 0,000013396 
B1 : integer := -14105;          -- = ((2^31)/1,995) * -0,000013103
B2 : integer := 14419;         -- = ((2^31)/1,995) * 0,000013396 
A1 : integer := -2147268361;        -- = ((2^31)/1,995) * -1,9948 
A2 : integer := 1070803162          -- = ((2^31)/1,995) * 0,99477

Die Koeffizienten sind in der Matlab/Octave Datei erzeugbar und können 
aus der Variable SOS herausgelesen werden.

[Martin O. (ossi-2)]

Hier mal ein Filter (in LUA) das gut für Debugging Zwecke geeignet ist, 
weil bei den ersten 20 Werten noch kein Fehler durch die Rundung 
passiert.
Alle Werte kann man leicht durch Handrechnung nachvollziehen.

scale=65536
B0=1.0*scale
B1=1.0*scale
B2=-0.5*scale
A1=-1.0*scale
A2=0.5*scale

x1=0
x2=0
y1=0
y2=0

for k=0,20 do
  if(k==2) then IIRinput=6400 else IIRinput=0 end
  x0=IIRinput

  yNew=(B0*x0+B1*x1+B2*x2-A1*y1-A2*y2)/scale
  x2=x1 ; x1=x0 ; y2=y1 ; y1=yNew

  print(k,IIRinput,yNew)
  end


k     input   output
0       0       0
1       0       0
2      6400     6400
3       0       12800
4       0       6400
5       0       0
6       0       -3200
7       0       -3200
8       0       -1600
9       0        0
10      0       800
11      0       800
12      0       400
13      0       0
14      0       -200
15      0       -200
16      0       -100
17      0       0
18      0       50
19      0       50
20      0       25

[Edi M. (Gast)]

Wozu benötigt, das das / 1,995? Um den Faktor 2 und einen Überlauf zu 
meiden?

[Marcel D. (diablokiller999)]

Meister E. schrieb:
> Wozu benötigt, das das / 1,995? Um den Faktor 2 und einen Überlauf zu
> meiden?

Um den Wertebereich auf 32 Bit (31 wegen signed) abzubilden und nicht 
mit Floats auf dem FPGA rechnen zu müssen. Ich normalisiere mehr oder 
weniger, heißt ich setze die 1,95 (nahe am höchsten Wert) als 2^31 fest 
und bekomme so einen Multiplikator, den ich auf alle Zahlen anwende.

[Martin O. (ossi-2)]

Du skalierst die Koeffizienten durch Multiplikation mit (2^31)/1,995.
Dann dividierst Du bei nA2(nA2'left downto INPUT_WIDTH) aber durch
2^32. Das passt nicht zueinander. Ich vermute aber, dass eventuell noch 
andere Fehler im Code sind.

Du könntest den Code eventull mal nach C oder JAVA portieren und dort 
testen, das finde ich immer viel einfacher als HDL Code zu debuggen.

[Detlef _. (detlef_a)]

Das Filter ist instabil. So funzt das nich.

>> roots([1 -1.9948 0.99477 ])

ans =

    1.0035
    0.9913


Cheers
Detlef

[Gerald M. (gerald_m17)]

Wenn du mit matlab schon deinen Filter erstellt hast, kannst du ihn auch 
problemlos mit dem HDL Coder inkl. Testbench als VHDL Datei speichern.

Wenn es dir um die Übung an sich geht, dann immer weiter.

[Gustl B. (-gb-)]

Mit pyFDA https://github.com/chipmuenk/pyFDA/blob/master/README.md kann 
man sich auch sehr viel schön einfach berechnen lassen. Vor allem auch 
die Koeffizienten in verschiedenen Wertebereichen.

[Marcel D. (diablokiller999)]

Gerald M. schrieb:
> Wenn du mit matlab schon deinen Filter erstellt hast, kannst du ihn auch
> problemlos mit dem HDL Coder inkl. Testbench als VHDL Datei speichern.
Ist in Octave geschrieben, glaube nicht, dass es dafür ein Plugin für 
sowas gibt?

Detlef _. schrieb:
> Das Filter ist instabil. So funzt das nich.

Laut PN-Plan liegen meine Pole innerhalb des Einheitskreises (siehe 
Anhang), soweit ich mich ans Studium erinnere, ist das Filter doch dann 
stabil oder nicht?

[Martin O. (ossi-2)]

Ich erhalte das gleiche Ergebnis wie Detlef _a für Deine Koeffizienen,
nämlich dass das Filter instabil ist. Aber vermutlich ist das nicht das 
einzige Problem.

[Detlef _. (detlef_a)]

Marcel D. schrieb:
> Gerald M. schrieb:
>> Wenn du mit matlab schon deinen Filter erstellt hast, kannst du ihn auch
>> problemlos mit dem HDL Coder inkl. Testbench als VHDL Datei speichern.
> Ist in Octave geschrieben, glaube nicht, dass es dafür ein Plugin für
> sowas gibt?
>
> Detlef _. schrieb:
>> Das Filter ist instabil. So funzt das nich.
>
> Laut PN-Plan liegen meine Pole innerhalb des Einheitskreises (siehe
> Anhang),

Nein

soweit ich mich ans Studium erinnere, ist das Filter doch dann
> stabil oder nicht?

Ja.

z^2 -1.9948*z + 0.99477 = 0

Davon die Nullstellen berechnen. Das geht ohne Programm mit der 
pq-Formel.

math rulez!
schönen Tach
Cheers
Detlef

[Marcel D. (diablokiller999)]

Detlef _. schrieb:
> z^2 -1.9948*z + 0.99477 = 0
>
> Davon die Nullstellen berechnen. Das geht ohne Programm mit der
> pq-Formel.
>

In meinem PN-Plan wird mit genaueren Werten gerechnet als es die SOS 
Struktur widerspiegelt, daher wohl die Instabilität.

>>> disp (sos);
>> sos =
>>  1.3396e-005  -1.3103e-005  1.3396e-005  1.0000e+000  -1.9948e+000  9.9477e-001
>> a
>> a =
>>  1.00000  -1.99476   0.99477

Aber dennoch bekomme ich mein Filter nicht zum laufen, ist meine VHDL 
Implementierung prinzipbedingt falsch? Ich glaube, der Hund liegt im 
shiften/cutten der Bits begraben -_-

[Detlef _. (detlef_a)]

>>>>>
In meinem PN-Plan wird mit genaueren Werten gerechnet als es die SOS
Struktur widerspiegelt, daher wohl die Instabilität.

Ja, durch die Normierung z.B. auf 32Bit wandern die Pole und das Ding 
kann instabil werden. Das ist das Hauptproblem bei 'scharfen' IIR deren 
Pole nahe am Einheitskreis liegen. Man braucht hohe Genauigkeit der 
Koeffizienten und damit zusammenhängend eine hohe Rechendynamik, 64Bit 
reicht da vllt. nicht.

Du trennst die Probleme des Entwurfs von den Problemen der 
Implementierung: Erst in Octave oder Scilab oder Matlab das Filter zum 
Fliegen bringen mit allen Koeffizientenrundungen und 
Bitbreitenbeschränkungen, dann die Migration aufs FPGA.

Schöne Aufgabe, viel Spass.
Cheers
Detlef

[Martin O. (ossi-2)]

Mit den Filterkoeffizienten von meinem oben gezeigten IIR Filter ist das
Debugging wegen der einfachen Zahlen einfacher. Manchmal hilft es auch, 
wenn
man zuerst in C das ganze so laufen lässt, wie es später in VHDL laufen 
soll.
In C ist das Debuggen viel einfacher.

[Marcel D. (diablokiller999)]

Detlef _. schrieb:
> Du trennst die Probleme des Entwurfs von den Problemen der
> Implementierung: Erst in Octave oder Scilab oder Matlab das Filter zum
> Fliegen bringen mit allen Koeffizientenrundungen und
> Bitbreitenbeschränkungen, dann die Migration aufs FPGA.
>

Da ich ja nicht an ein gewisses Filterverhalten gebunden bin, hab ich 
mir mal ein stabiles Butterworth-Filter mit 2 Koeffizienten gebastelt:

SOS : 0.097631   0.195262   0.097631   1.000000  -0.942809   0.333333
roots([1.000000  -0.942809   0.333333])
=> 0.47140 + 0.33333i
=> 0.47140 - 0.33333i

Zudem habe ich mich doch mal an eine Implementierung mit asynchronem 
Verhalten getraut (nach transponierter DF2). Mit Filterkoeffizienten aus 
einem Buch sieht das ganze schon nach einem normalen Filterverhalten aus 
- zumindest übersteuert es nicht wie blöde. Meine errechneten 
Koeffizienten bringen das Teil aber anscheinend wieder durcheinander. 
Irgendwo ist noch der Wurm drin...

[Martin O. (ossi-2)]

Ich finde FPGA-Code, der komplizierte Rechnungen macht, immer schwer zu 
debuggen. Ich mache dann oft erstmal eine "Simulation" des Codes in 
einer passenden Sprache. Das hab ich mit Deinem VHDL Code gemacht. Als 
LUA Code sieht das ganze dann so aus

scale=65536

cB0= 1.0*scale
cB1=-2.0*scale
cB2= 3.0*scale
cA1=-1.0*scale
cA2= 0.5*scale

nZ1=0
nZ2=0

for k=0,10 do
  if(k==1) then IIRinput=800 else IIRinput=0 end
  nX=IIRinput
  print("nX=",nX)
--mul
  nA1 = nZ1 * cA1;
  nA2 = nZ2 * cA2;
  nB1 = nZ1 * cB1;
  nB2 = nZ2 * cB2;
  print("nA1=",nA1/scale," nA2=",nA2/scale)
--s1
  nSUMA1A2 = nA1/scale + nA2/scale;
  nSUMB1B2 = nB1/scale + nB2/scale;
  print("nSUMA1A2=",nSUMA1A2," nSUMB1B2=",nSUMB1B2 )
--s2
  nSUMX = nX -  nSUMA1A2;
  print("nSUMX=",nSUMX)
--s3;
  nZ1help = nSUMX;
  nZ2help = nZ1;
  nZ1=nZ1help
  nZ2=nZ2help
  nB0 = nSUMX * cB0;
--s4;
  nSUMXB0 = nB0/scale + nSUMB1B2;
  print("nSUMXB0 ==============================================================",nSUMXB0 ) ;
--convert
   --nY          <= std_logic_vector(nSUMXB0(nSUMXB0'left downto nSUMXB0'left-INPUT_WIDTH+1));
--finished
  --oIIR_TX <= nY(nY'left downto (nY'left - OUTPUT_WIDTH+1));
  end

Der Code  ist ziemlich nahe am Original. Die Bitshifts hab ich durch 
Divisionen durch den Skalierungsfaktor ersetzt. Ich rechne ein einfaches 
Filter, das man besser von Hand nachrechnen kann. Damit habe ich Deinen 
Code analysiert und folgenden Fehler gefunden:

Die Zeile

 nSUMA1A2    <= nA1(nA1'left downto INPUT_WIDTH) - nA2(nA2'left downto 
INPUT_WIDTH);

muss ersetzt werden durch:

 nSUMA1A2    <= nA1(nA1'left downto INPUT_WIDTH) + nA2(nA2'left downto 
INPUT_WIDTH);


Dann ist, glaube ich, die Skalierung des Ausgangssignals unnötig (vgl. 
LUA Code). Weiter musst Du aufpassen (hatte ich oben schon gesagt): Dein
Code verwendet nach den Multiplikationen eine Division durch 2 ^ 
Input_Width, damit musst Du auch deine Koeffzienten skalieren (nicht mit 
1.95....)

Zur Simulation Deines Filters kannst Du ja den LUA Code nehmen.

[Marcel D. (diablokiller999)]

Martin O. schrieb:
> Weiter musst Du aufpassen (hatte ich oben schon gesagt): Dein
> Code verwendet nach den Multiplikationen eine Division durch 2 ^
> Input_Width, damit musst Du auch deine Koeffzienten skalieren (nicht mit
> 1.95....)
>
> Zur Simulation Deines Filters kannst Du ja den LUA Code nehmen.

Wow, vielen Dank für die Mühe!
Es hat den Anschein, als habe ich die Umrechnung der Koeffizienten noch 
nicht ganz verstanden und genau dort vermute ich gerade den Fehler, 
nachdem meine weiteren Versuche einen Filter zu basteln ähnlich 
scheitern wie mit diesem Code hier. Wäre super wenn mich da jemand 
erleuchten könnte, wie ich mein 48Bit Eingangssignal mit den 
Matlab-Koeffizienten verrechnen muss :)

[Edi M. (Gast)]

Ja, das glaube ich auch. Mit den Koeffizienten ist definitiv was faul!

[Martin O. (ossi-2)]

Ein paar Bemerkungen zur Skalierung etc.
Zum Üben kann man das folgende Filter nehmen:

scale=65536
gain=0.075
B0=1.0*gain*scale
B1=2.0*gain*scale
B2=1.0*gain*scale
A1=-1.4*scale
A2= 0.7*scale

x1=0
x2=0
y1=0
y2=0

for k=0,28 do
  if(k>=2) then IIRinput=5000  else IIRinput=0 end
  x0=IIRinput

  yNew=(B0*x0+B1*x1+B2*x2-A1*y1-A2*y2)/scale
  x2=x1 ; x1=x0 ; y2=y1 ; y1=math.floor(yNew)

  print(k,IIRinput,y1)
  end

k     Input     y1
00  0
1  0  0
2  5000  375
3  5000  1650
4  5000  3547
5  5000  5310
6  5000  6451
7  5000  6814
8  5000  6523
9  5000  5862
10  5000  5140
11  5000  4592
12  5000  4330
13  5000  4347
14  5000  4554
15  5000  4832
16  5000  5077
17  5000  5225

Es ist ein Tiefpassfilter mit näherungsweise Tschebyscheff 
Charakteristik.

Gestartet wurde das Design mit den Filterkoeffizienten
b0=1 ; b1=2 ; b2=-1 ; a1=-1.4 ; a2=0.7
Diese ergeben sich aus den Pol- und Nullstellen. Dieses Filter hat für 
Frequenz 0 allerdings eine Verstärkung von A=13.3333=1/0.075. Um ein 
Filter zu erhalten das bei f=0 die Verstärkung 1 hat, multipliziert man 
die b-Koeffizienten mit 0.075=1/A.

Nun sollen die Koeffizienten skaliert werden. Es werden alle 
Koeffizienten mit S=65536=2^16 multipliziert. Hinter den Multilizierern 
wird dann durch S dividiert, indem die Resultate um 16 nach rechts 
schiebt. Soweit ist es noch simpel. Der kompliziertere Teil kommt 
danach:

Jetzt muss man festlegen, mit welchen Wortbreiten man addiert und 
multipliziert. Die Bitbreiten muss man so bestimmen, dass kein Überlauf 
stattfinden kann. Die Bitbreite des Eingangssignals ist meist bekannt.
Die Bitbreite des Ausgangssignals meist auch. Schwierig sind
die Bitbreiten der internen Größen. Da IIR Filter resonantes Verhalten 
zeigen können, können diese internen Signale wesentlich größer werden, 
als die Eingangs- uind Ausgangssignale. Die "interne Verstärkung" kann 
man,  glaube ich, anhand der Polstellen näherungsweise bestimmen. Ich 
nehme aber oft einfach 4 Bits mehr als unbedingt nötig, und mache dann 
Tests mit verschiedensten Eingangssignalen. Vielleicht kennt hier jemand 
ein besseres Vorgehen.

[Marcel D. (diablokiller999)]

Danke für die ausführliche Erklärung :)
Ich habe anscheinend einen laufenden IIR auf dem FPGA implementiert, 
werde ihn heute noch versuchen zu testen aber die Simulation schaut 
vernünftig aus. Nur klappt das mit der Verstärkung von 1 nicht genau, 
was aber wohl an Quantisierungseffekten liegen dürfte.
Meine Ergebnisse und den Sourcecode/Testbench lasse ich hier mal für 
andere stehen, eventuell hilft es ja jemand anderem Ratlosen weiter.

[Martin O. (ossi-2)]

Was für ein Filter soll das denn sein? Woher stammt das Design?
Wie lauten die Filterkoeffizienten unskaliert?

[Marcel D. (diablokiller999)]

Ein einfacher 2-Stage Butterworth Filter zum testen des SOS-Designs. Ich 
habe ihn im FPGA hinter ein 4 Stage CIC Filter verbaut. Die 
Koeffizienten für das SOS-Design sind

> b0 = 0.097631
> b1 = 0.195262
> b2 = 0.097631
> a0 = 1.000000
> a1 = -0.942809
> a2 = 0.333333

Die Koeffizienten habe ich mit (2^31/0,943)*Wert skaliert, der Gain 
müsste bei 1,00000085270047 liegen, wenn ich mich nicht verrechnet habe. 
Simulation und FPGA-Test zeigen in etwa das gleiche Bild, somit sollte 
es geklappt haben.

[Martin O. (ossi-2)]

Ich glaube in Deiner aktuellen Version steckt wieder ein Fehler
Statt
nSUMX <= nX + nSUMA1A2(67 downto 0);
muss es heissen
nSUMX <= nX - nSUMA1A2(67 downto 0);
damit die Berechnung mit dem Blockplan aus dem allerersten Posting 
übereinstimmt.

Weiter kannst Du nicht die Koeffizienten mit (2^31/0,943)
skalieren, und die Multiplikationsergebnisse durch 2^31 dividieren.
Die beiden Skalenfaktoren müssen gleich sein, falls nicht, änderst
Du die Filtercharakteristik.

Du solltest mal ein Filter laufen lassen, bei dem Du jedes Bit
genau nachvollziehen kannst. Ich habe ja Beispiele dazu gegeben.
Alles andere ist ein Stochern im Nebel......

[Detlef _. (detlef_a)]

Hi, deine Eingangsdaten sind 64 Bit breit, scale macht 17Bit, gain 
verliert wieder 4Bit, 4 Additionen je 1 Bit, engineering margin 4Bit, 
macht 85 Bit nötige Wortbreite.

So breit rechnest Du nicht?!

Cheers
Detlef
? Jetzt ist der Beitrag weg, auf den ich mich bezog. Schade.

[Marcel D. (diablokiller999)]

Martin O. schrieb:
> Statt
> nSUMX <= nX + nSUMA1A2(67 downto 0);
> muss es heissen
> nSUMX <= nX - nSUMA1A2(67 downto 0);

Ich glaube das + ist schon richtig, da der Wert von A1 negativ ist.

Ich habe meine Implementierung mal auf 16Bit Koeffizienten minimiert und 
gebe Werte aus dem LUA-Skript ein. Dennoch verhalten sie sich nicht 
gleich, wenn ich das + durch ein - ersetze wird alles nur noch 
schlimmer. Hoffe die Kommentare geben Aufschluss über die 
Implementierung, ich multipliziere meine Koeffizienten mit 2^16 und 
right shifte um 2^16 nach der Multiplikation.

scale=65535
gB0=1
gB1=2*gB0
gB2=1
gA1=-1.4
gA2=0.7
gain = 1/((gB0+gB1+gB2)/(1+gA1+gA2))
B0=gB0*gain*scale
B1=gB1*gain*scale
B2=gB2*gain*scale
A1=gA1*scale
A2= gA2*scale
print(B0,B1,B2,A1,A2,gain)
x1=0
x2=0
y1=0
y2=0
-- Input 64Bit signed
for k=0,60 do
  if(k<=20 or k>=40) then IIRinput=9223372036854775807  else IIRinput=0 end
  x0=IIRinput
  
  yNew=(B0*x0+B1*x1+B2*x2-A1*y1-A2*y2)/scale
  x2=x1 ; x1=x0 ; y2=y1 ; y1=math.floor(yNew)

  print(k,IIRinput,y1)
  end

Detlef _. schrieb:

> ? Jetzt ist der Beitrag weg, auf den ich mich bezog. Schade.

Hab was editiert, aber man kann hier ja leider keine Anhänge tauschen ^^

Detlef _. schrieb:
> Hi, deine Eingangsdaten sind 64 Bit breit, scale macht 17Bit, gain
> verliert wieder 4Bit, 4 Additionen je 1 Bit, engineering margin 4Bit,
> macht 85 Bit nötige Wortbreite.
>
> So breit rechnest Du nicht?!

signal nB0      : signed(83 downto 0) := (others => '0'); --  64 + 4 Guard + 16 Bit
signal nB1      : signed(83 downto 0) := (others => '0'); --  64 + 4 Guard + 16 Bit
signal nB2      : signed(83 downto 0) := (others => '0'); --  64 + 4 Guard + 16 Bit
signal nA1      : signed(84 downto 0) := (others => '0'); --  64 + 4 Guard + 17 Bit
signal nA2      : signed(83 downto 0) := (others => '0'); --  64 + 4 Guard + 16 Bit

Sind doch 85/84 Bit für die Multiplikationen. Die Addierer können ja 
kleiner sein, weil ich vorher eh shifte?

[Detlef _. (detlef_a)]

Ok, hatte ich nicht wahrgenommen. Dein oIIR_Tx sieht erwartungskonform 
aus,  6.28/angle(roots([1 -1.4 0.7])) ~ 11, dh. dein Filter schwingt 
gedämpft mit Periodendauer 11 Abtastwerte. Das passt doch, die 
Vorzeichen sind richtig.

>>>Dennoch verhalten sie sich nicht gleich,

Was verhält sich nicht gleich. Sieht alles gut aus. Was geht jetzt 
nicht?

Cheers
Detlef

[Martin O. (ossi-2)]

Ich wette um 10 cent dass hier das Minuszeichen hinmus

nSUMX <= nX - nSUMA1A2(67 downto 0);

genauso wie hier vor den A Produkten jeweils Minuszeichen sind....

yNew=(B0*x0+B1*x1+B2*x2-A1*y1-A2*y2)/scale

[Marcel D. (diablokiller999)]

Detlef _. schrieb:
> Ok, hatte ich nicht wahrgenommen. Dein oIIR_Tx sieht erwartungskonform
> aus,  6.28/angle(roots([1 -1.4 0.7])) ~ 11, dh. dein Filter schwingt
> gedämpft mit Periodendauer 11 Abtastwerte. Das passt doch, die
> Vorzeichen sind richtig.
Kannst du das genauer erklären? :D

> Was verhält sich nicht gleich. Sieht alles gut aus. Was geht jetzt
> nicht?

Müsste ich durch meine Skalierung mit dem Verstärkungsfaktor bei f0 
nicht das selbe Ausgangs- wie Eingangssignal haben? Zudem sind die 
Ausgangswerte nicht identisch mit dem Lua-Skript...

Martin O. schrieb:
> Ich wette um 10 cent dass hier das Minuszeichen hinmus
>
> nSUMX <= nX - nSUMA1A2(67 downto 0);
>

Dann pendelt sich der Wert jedenfalls nicht wieder bei 0 ein in der 
Simulation (Siehe Bild). Hab meinen Testbench ja mit beigelegt, könnt 
simulieren wenn ihr wollt :)

[Detlef _. (detlef_a)]

In deinem letzten post sieht das Ausgangssignal aber anders aus als in 
dem von 12:03, auf das ich mich bezogen hatte.

Marcel D. schrieb:
> Detlef _. schrieb:
>> Ok, hatte ich nicht wahrgenommen. Dein oIIR_Tx sieht erwartungskonform
>> aus,  6.28/angle(roots([1 -1.4 0.7])) ~ 11, dh. dein Filter schwingt
>> gedämpft mit Periodendauer 11 Abtastwerte. Das passt doch, die
>> Vorzeichen sind richtig.
> Kannst du das genauer erklären? :D
>

Die Pole des Filters sind bei roots([1 -1.4 0.7]) ~0.7+-j0.46. Die 
Zahlen haben den Winkel +-0.57rad, 2*pi rad ist die Abtastfrequenz, das 
ist ~1/11 davon.

>> Was verhält sich nicht gleich. Sieht alles gut aus. Was geht jetzt
>> nicht?
>
> Müsste ich durch meine Skalierung mit dem Verstärkungsfaktor bei f0
> nicht das selbe Ausgangs- wie Eingangssignal haben?

Du gehst rein mit einem Sprung, dann siehst du die Sprungantwort des 
Filters, und die sieht qualitativ richtig aus, auf die Amplitude habe 
ich nicht geschaut.

Das Vorzeichen stimmt, die Wette halte ich, siehe z.B. slide 23: 
http://public.beuth-hochschule.de/~purat/lehre/dsv2/slides/Realisierung.pdf

Cheers
Detlef

[Marcel D. (diablokiller999)]

Detlef _. schrieb:
> In deinem letzten post sieht das Ausgangssignal aber anders aus als in
> dem von 12:03, auf das ich mich bezogen hatte.
>
Das war auch der Versuch mit Ossis Minus :)

> Du gehst rein mit einem Sprung, dann siehst du die Sprungantwort des
> Filters, und die sieht qualitativ richtig aus, auf die Amplitude habe
> ich nicht geschaut.

Genau das macht mir halt noch Bauchschmerzen, von dem was ich hier 
gelernt habe, müsste die Amplitude meinem Eingangssprung entsprechen, 
wenn ich bei der Skalierung alles richtig gemacht habe. Die Summe aus 
B1+B2+B3 dividiert durch A1+A2+A3 ergibt für mich 13,33 und mein Faktor 
für B1 bis B3 also 0,0075 (siehe Lua skript).

> Das Vorzeichen stimmt, die Wette halte ich, siehe z.B. slide 23:
> http://public.beuth-hochschule.de/~purat/lehre/dsv2/slides/Realisierung.pdf
>
Nice, was zum schmöckern :)

[Martin O. (ossi-2)]

Das Vorzeichen stimmt, die Wette halte ich, siehe z.B. slide 23:
http://public.beuth-hochschule.de/~purat/lehre/dsv...

Da haben die A-Koeffizienten genau ein Minuszeichen,
genau wie ichs meine.

[Weltbester FPGA-Pongo (Gast)]

Was sagt denn die Simulation? Vielleicht hast Du es in VHDL nur nicht 
richtig hingeschrieben, eine Klammer vergessen, einen Datentypen falsch 
bezeichnet oder eine Konstante negiert.

[Marcel D. (diablokiller999)]

Weltbester FPGA-Pongo schrieb im Beitrag #5119151:
> Was sagt denn die Simulation? Vielleicht hast Du es in VHDL nur nicht
> richtig hingeschrieben, eine Klammer vergessen, einen Datentypen falsch
> bezeichnet oder eine Konstante negiert.

Sollte soweit alles stimmen, der Code und die Simulation habe ich ja 
inklusive Screenshots gepostet :)

[Weltbester FPGA-Pongo (Gast)]

Wo ist dann das Problem?

[Detlef _. (detlef_a)]

Martin O. schrieb:
> Das Vorzeichen stimmt, die Wette halte ich, siehe z.B. slide 23:
> http://public.beuth-hochschule.de/~purat/lehre/dsv...
>
> Da haben die A-Koeffizienten genau ein Minuszeichen,
> genau wie ichs meine.

>>>>>>>>yNew=(B0*x0+B1*x1+B2*x2-A1*y1-A2*y2)/scale

hamse doch.
Cheers
Detlef

[Weltbester FPGA-Pongo (Gast)]

Ok, wobei man an der Schreibweise noch feilen könnte.

[Marcel D. (diablokiller999)]

Weltbester FPGA-Pongo schrieb im Beitrag #5119187:
> Wo ist dann das Problem?

Das Problem ist, dass ich kein konsistentes Verhalten bei meiner 
Implementierung finden kann. Je nachdem, ob ich ein Minus

nSUMA1A2 <= nA1(84 downto 16) - nA2(83 downto 16);

oder Plus

nSUMA1A2 <= nA1(84 downto 16) + nA2(83 downto 16);

an dieser Stelle implementiere, bekomme ich ein unterschiedliches 
Verhalten von meinen Filterkoeffizienten
> b0 = 0.097626
> b1 = 0.195251
> b2 = 0.097626
> a1 = -0.942810
> a2 = 0.333328
> Gain = 1.00003841199684
> siehe Own Coeff Minus/Plus

und den Koeffizienten von Ossi
> b0 = 1
> b1 = 2
> b2 = 1
> a1 = -1.4
> a2 = 0.7
> Gain = 0.0075
> siehe Random Coeff Minus/Plus
Alle B-Koeffizienten sind dabei mit dem Gain multipliziert.

Wenn ich mich richtig an meine Vorlesung erinnere, müssen alle 
rückgekoppelten Signale subtrahiert werden, also wäre doch eher folgende 
Implementierung richtig

nSUMA1A2 <= nA1(84 downto 16) - nA2(83 downto 16);  -- divide by 2^16
nSUMB1B2 <= nB1(83 downto 16) + nB2(83 downto 16);  -- divide by 2^16
nSUMXB0 <= nB0(83 downto 16) + nSUMB1B2;            -- divide by 2^16
nSUMX <= nX - nSUMA1A2(67 downto 0);

Dabei sehen die Kurven in "Own Coeff Minus" und "Random Coeff Plus" für 
mich qualitativ richtig aus, zumindest pendeln sie sich um einen Wert 
ein und finden auch wieder zur 0 zurück. Dabei muss aber eine 
Implementierung falsch sein und ich blicke langsam echt nicht mehr 
durch.

Detlef _. schrieb:
>> Da haben die A-Koeffizienten genau ein Minuszeichen,
>> genau wie ichs meine.
>
>>>>>>>>>yNew=(B0*x0+B1*x1+B2*x2-A1*y1-A2*y2)/scale
>
> hamse doch.
> Cheers
> Detlef

Mein A2 ist aber kein negativer Koeffizient ^^

[Martin O. (ossi-2)]

Damit es mit den Vorzeichen keine Verwirrung gibt gehe ich üblicherweise
in meinen Programmen wie folgt vor.

Die Übertragungsfunktion des biquads ist

$$H(z)={ b_0 + b_1 z^{-1} + b_2 z^{-2} \over 1 + a_1 z^{-1} + a_2 z^{-2}}$$

Und diese Bezeichnungen gelten, egal ob die einzelnen Koeffizienten
positiv oder negativ sind. In meinen Programmen stehen dann die 
Variablen
B0,B1,... für genau diese Koeffizienten und sie können im einzelnen
wieder positiv und negativ sein. Mit diesen Variablen gilt dann
die Rekursion

 yNew=(B0*x0+B1*x1+B2*x2-A1*y1-A2*y2)

wobei vor den As jetzt Minuszeichen auftauchen. Das Filter programmiert
man dann ohne dass man weiss, ob die einzelnen Koeffizienten positiv 
oder
negativ sind.

nSUMA1A2 <= nA1(84 downto 16) - nA2(83 downto 16);  -- divide by 2^16
nSUMB1B2 <= nB1(83 downto 16) + nB2(83 downto 16);  -- divide by 2^16
nSUMXB0 <= nB0(83 downto 16) + nSUMB1B2;            -- divide by 2^16
nSUMX <= nX - nSUMA1A2(67 downto 0);

Dieser Block ist dann falsch, weil A1 und A2 mit unterschiedlichen
Vorzeichen verrechnet werden. Meiner Meinung nach kannst Du schreiben:

nSUMA1A2 <= nA1(84 downto 16) + nA2(83 downto 16);  -- divide by 2^16
nSUMB1B2 <= nB1(83 downto 16) + nB2(83 downto 16);  -- divide by 2^16
nSUMXB0 <= nB0(83 downto 16) + nSUMB1B2;            -- divide by 2^16
nSUMX <= nX - nSUMA1A2(67 downto 0);

oder

nSUMA1A2 <= -nA1(84 downto 16) - nA2(83 downto 16);  -- divide by 2^16
nSUMB1B2 <= nB1(83 downto 16) + nB2(83 downto 16);  -- divide by 2^16
nSUMXB0 <= nB0(83 downto 16) + nSUMB1B2;            -- divide by 2^16
nSUMX <= nX + nSUMA1A2(67 downto 0);

Alle Klarheiten beseitigt?

[Martin O. (ossi-2)]

Womit simuliert Ihr VHDL Code ? Gibt es kostenlose VHDL Simulatoren?

[Michael W. (Gast)]

Martin O. schrieb:
> Womit simuliert Ihr VHDL Code ? Gibt es kostenlose VHDL Simulatoren?
Der ModelSIM von Altera ist frei.

[J. S. (engineer)]

Ich habe jetzt nicht alles angesehen, aber dafür, dass Du mit bis zu 64 
Bits rechnest, ist das Ergebnis von der ModelSIM-Grafik her 
"bescheiden".

Sicher, dass das vom VHDL her passt? Z.B. sollte man darüber ...

>nSUMA1A2 <= nA1(99 downto 32) + nA2(99 downto 32);
... nochmal nachdenken und auch generell über die Themen "Rundung" und 
"Fehlerfortplanzung".

Simuliere Dir das mal in Real und dann in INT. Wenn es in Real läuft, 
ist es entweder ein Rundungs- oder Timingproblem, also eine Frage der 
Zuordnung von Ergebnissen.

[Niklas (Gast)]

Martin O. schrieb:
> Womit simuliert Ihr VHDL Code ? Gibt es kostenlose VHDL Simulatoren?

GHDL + gtkwave

Ist zwar nur mit 'Fachwissen' zu benutzen, aber kostenlos & open 
source...

[Marcel D. (diablokiller999)]

So, meine neueste Implementierung. Getaktet und funktionierend, 
inklusive 1er Gain und Überlaufschutz. Mein Problem war wohl, dass ich 
nach der Multiplikation direkt dividierte und somit ein so großer Teil 
der Auflösung flöten ging, dass augenscheinlich nur Murks heraus kam. 
Sobald ich anfing, die Koeffizienten nur um 2^7 zu dividieren, mit dem 
genaueren Signal zu rechnen und den Rest der 2^11 erst vor den 
Koeffizienten wieder zu dividieren, hat die Implementierung 
funktioniert. Ossi hatte auch recht was die Vorzeichen angeht, danke 
dafür nochmal!

Bis auf die Zeile

nX      <= signed(iIIR_RX) & "0000000000000000000000";

sollte die Implementierung für jedes Q-Format und jede Bitweite 
funktionieren. Jedenfalls scheint meine Simulation soweit zu klappen :)

Zur Simulation nutze ich die kostenlose ModelSIM Version. Etwas 
umständlich, aber das ist in der FPGA Welt eh fast alles. Kann Xilinx 
mit ISE/Vivado aber empfehlen, hat mir persönlich mehr Spaß gemacht als 
Altera :)

[Marcel D. (diablokiller999)]

Wenn man die States Idle und Mul editiert:

when idle =>
    oDone       <= '0';
    if(iNewValue = '1') then                
        state   <= mul;
        nX(INPUT_WIDTH-1 downto 0) <= signed(iIIR_RX);
        
    end if;
when mul =>                     -- Multiply signals for sums
    nA1         <= nZ1(nZ1'left downto QFORMAT) * cA1;
    nA2         <= nZ2(nZ2'left downto QFORMAT) * cA2;
    nB1         <= nZ1(nZ1'left downto QFORMAT) * cB1;
    nB2         <= nZ2(nZ2'left downto QFORMAT) * cB2;
    nX      <= shift_left(signed(nX),QFORMAT);
    state       <= s1;

Sollte das gesamte Design leicht editierbar sein. Bin schon gespannt, 
wie sich das Filter im FPGA nun verhält und ob man es kaskadieren kann. 
Dann gehts ans optimieren ^^

[Marcel D. (diablokiller999)]

Nochmal ein Vergleich Matlab vs. FPGA
Scheint nun endlich zu funktionieren, ein dickes Danke nochmal an alle!
Habe mal ein paar Fehler aus der VHDL Datei behoben und alles etwas 
verschlankt, hoffe das einige Leute damit was anfangen können! :)

[Detlef _. (detlef_a)]

Die zerissenen 'Kuppen' des IIR und des CIC schmecken nach overflow. Du 
hast die Bitbreiten ja ganz prima parametrisiert, zieht die nen bißchen 
hoch und schaue, ob der Effekt weg ist. Ist zwar sowieso egal, weil das 
ja im Stopband liegt, sieht aber aus wie 'übertünchter Pfusch', hätte 
ich kein gutes Gefühl bei.

Cheers
Detlef

[Marcel D. (diablokiller999)]

Detlef _. schrieb:
> Die zerissenen 'Kuppen' des IIR und des CIC schmecken nach overflow.

Eigentlich habe ich einen Overflow-Schutz implementiert:

if((nSUMXB0(nSUMXB0'left downto nSUMXB0'left-5) = "000000") OR (nSUMXB0(nSUMXB0'left downto nSUMXB0'left-5) = "111111")) then
    oIIR_TX     <= std_logic_vector(nSUMXB0(nSUMXB0'left-4 downto 22));
else
    if(nSUMXB0(nSUMXB0'left) = '1') then
        oIIR_TX  <= (others=>'0');
        oIIR_TX(oIIR_TX'left) <= '1';
    else
        oIIR_TX  <= (others=>'1');
        oIIR_TX(oIIR_TX'left) <= '0';
    end if;
end if;

Die "Kuppen" scheinen auch eher vom vorgeschalteten CIC Filter zu 
kommen, da dieser genau dieses Muster aufweist. Die Kurve vom IIR ist 
nur errechnet aus dem Ausgangssignal - Signal nach CIC.

[Martin O. (ossi-2)]

Ich hätte einen Verbesserungsvorschlag:

Momentan benutzt Du den Wert QFORMAT für zwei (verschiedene) Zwecke:
a) Anzahl der Bits zur Koeffizientenspeicherung
b) Anzahl der Bits zun Schieben nach der Multiplikation
Diese beiden Werte müssen aber nicht unbedingt gleich sein.

Könntest Du eine Version erstellen, bei der man für die beiden Zwecke
getrennte Parameter (z.B. QSIZE und QSHIFT) hat. Leider reichen
meine VHDL Kenntnisse nicht aus, das selbst zu machen.

Mit der Modifikation kann man auch Koeffizienten B0,B1,B2 benutzen,
die grösser als 1 sind.

[Marcel D. (diablokiller999)]

Martin O. schrieb:
> Könntest Du eine Version erstellen, bei der man für die beiden Zwecke
> getrennte Parameter (z.B. QSIZE und QSHIFT) hat. Leider reichen
> meine VHDL Kenntnisse nicht aus, das selbst zu machen.
>
> Mit der Modifikation kann man auch Koeffizienten B0,B1,B2 benutzen,
> die grösser als 1 sind.

Der Shift muss identisch zum Q-Format sein, weil das Inputsignal und die 
Multiplikationsdaten die selbe Zehnerpotenz haben müssen,sonst haut das 
mit dem dividieren bei

nZ1         <= nSUMX(nSUMX'left downto QFORMAT);

nicht hin.

Was du möchtest ist das, was ich für nA1 gemacht habe. Soweit ich das 
sehe, musst du dafür "einfach" nur den Bitbereich für die anderen 
Multiplikationen und Konstanten erweitern.

constant cA1 : signed(QFORMAT+1 downto 0)  := to_signed(A1,QFORMAT+2);

und

signal nA1      : signed(INPUT_WIDTH-1+GUARDBITS+QFORMAT+2 downto 0);

sollten dafür die Blaupause sein.

Du müsstest also die ganze Deklaration wie folgt ändern:

constant cA1 : signed(QFORMAT+1 downto 0)  := to_signed(A1,QFORMAT+2);-- A1
constant cA2 : signed(QFORMAT+1 downto 0)  := to_signed(A2,QFORMAT+2);-- A2
constant cB0 : signed(QFORMAT+1 downto 0)  := to_signed(B0,QFORMAT+2);-- B1
constant cB1 : signed(QFORMAT+1 downto 0)  := to_signed(B1,QFORMAT+2);-- B1
constant cB2 : signed(QFORMAT+1 downto 0)  := to_signed(B2,QFORMAT+2);-- B1

signal nSUMX    : signed(INPUT_WIDTH-1+GUARDBITS+QFORMAT downto 0);             -- Sum of shifted Input and recursive sum
signal nSUMA1A2    : signed(INPUT_WIDTH-1+GUARDBITS+QFORMAT+1 downto 0);        -- Sum of recursive multiplication
signal nSUMB1B2    : signed(INPUT_WIDTH-1+GUARDBITS+QFORMAT+1 downto 0);          -- Sum of forward multiplication
signal nSUMXB0    : signed(INPUT_WIDTH-1+GUARDBITS+QFORMAT+1 downto 0);           -- Sum before Output
signal nB0      : signed(INPUT_WIDTH-1+GUARDBITS+QFORMAT+2 downto 0);           -- Multiplication of B0
signal nB1      : signed(INPUT_WIDTH-1+GUARDBITS+QFORMAT+2 downto 0);           -- Multiplication of B1
signal nB2      : signed(INPUT_WIDTH-1+GUARDBITS+QFORMAT+2 downto 0);           -- Multiplication of B2
signal nA1      : signed(INPUT_WIDTH-1+GUARDBITS+QFORMAT+2 downto 0);           -- Multiplication of A1
signal nA2      : signed(INPUT_WIDTH-1+GUARDBITS+QFORMAT+2 downto 0);           -- Multiplication of A2

und den Abgriff der Daten um ein "zusätzliches" Guardbit erweitern

                    
if((nSUMXB0(nSUMXB0'left downto nSUMXB0'left-(GUARDBITS+2)) = "0000000") OR (nSUMXB0(nSUMXB0'left downto nSUMXB0'left-(GUARDBITS+2)) = "1111111")) then
                        oIIR_TX     <= std_logic_vector(nSUMXB0(nSUMXB0'left-GUARDBITS-1 downto QFORMAT));

Meine Simulation funktioniert damit noch einwandfrei, ob das nun mit 
größeren Werten als 1 funktioniert, konnte ich leider noch nicht testen.

[Marcel D. (diablokiller999)]

Nun habe ich versucht, beide Filter hintereinander zu kaskadieren, 
leider ohne Erfolg. Meine MATLAB-Simulation scheint dabei stabil, auch 
die einzelnen A's der beiden Filter sind laut PN-Plan stabil. Was habe 
ich übersehen?

Meine Koeffizienten für den Filter 4. Ordnung laut tf2sos sind:
SOS1:
a1= -1,67E+00
a2 = 7,05E-01
b0 = 1,33E-04
b1 = 1,98E-04
b2 = 1,33E-04

SOS2:
a1 = -1,83E+00
a2 = 8,67E-01
b0 = 1
b1 = 1,58E-01
b2 = 1

B-Gain = 0,029321379760263

[Gustl B. (-gb-)]

Jo, die Koeffizienten sehen gut aus, funktioniert.

Hab das mal in Python geschrieben und werde das nachher mal in VHDL 
übersetzen. Für die Simulation bietet sich ein Sweep an der also schön 
linear die Frequenz vom Sinus als Eingangssignal erhöht, kann man auch 
prima in VHDL schreiben und simulieren.

[Martin O. (ossi-2)]

So, ich habe jetzt ne Version erstellt, bei der man die
Parameter:
a) QSIZE: Anzahl der Bits zur Koeffizientenspeicherung
b) QSHIFT: Anzahl der Bits zun Schieben nach der Multiplikation
getrennt einstellen kann. Dann habe ich die testbench so modifiziert,
dass Testausgaben erzeugt werden, mit denen man die Werte
innerhelb des Filters nachkontrollieren kann.
Ausgabe sieht wie folgt aus:

# nValue= 0     nOUT= 0     help 1..4 /65536= 0         0      0       0 at time 100 ns
# nValue= 0     nOUT= 0     help 1..4 /65536= 0         0      0       0 at time 200 ns
# nValue= 1000  nOUT= 0     help 1..4 /65536= 0         0      0       0 at time 300 ns
# nValue= 0     nOUT= 1000  help 1..4 /65536= 0         0      0       0 at time 400 ns
# nValue= 0     nOUT= 3000  help 1..4 /65536= -1000     0   2000       0 at time 500 ns
# nValue= 0     nOUT= 5500  help 1..4 /65536= -1000   500   2000    3000 at time 600   ns
# nValue= 0     nOUT= 4000  help 1..4 /65536=  -500   500   1000    3000 at time 700 ns

[Gustl B. (-gb-)]

Hm, wie würden das die Profis in VHDL machen, mit integer, signed, oder 
lieber alles mit real?

[J. S. (engineer)]

"signed", für real gibt es da keine Notwendigkeit.
"überlaufschutz" braucht es auch keinen, wenn die Dynamik stimmt
da muss man eben passend zwischen skalieren

[Gustl B. (-gb-)]

Wie macht man das dann mit den Koeffizienten?

#SOS1:
IIR_0_a1= -1.67E+00
IIR_0_a2 = 7.05E-01
IIR_0_b0 = 1.33E-04
IIR_0_b1 = 1.98E-04
IIR_0_b2 = 1.33E-04

#SOS2:
IIR_1_a1 = -1.83E+00
IIR_1_a2 = 8.67E-01
IIR_1_b0 = 1.0
IIR_1_b1 = 1.58E-01
IIR_1_b2 = 1.0

B_Gain = 0.029321379760263

Also wenn ich das mit signed(15 downto 0) mache, kann ich das einfach so 
skalieren, dass -1.67E+00 dann -2**15 ist? Oder eben +2=2**15-1 und 
-2=-2**15.
Oder muss ich bei Koefizienten <1 tatsächlich eine Division einbauen?

[J. S. (engineer)]

Das hängt jetzt ein wenig an den Gleichungen und dem Wert der 
Koeffizienten - ganz optimal, also schnell UND klein wird es wenn man 
sich für jeden Schritt die optimale Breite des Vektors rechnet, unter 
Berücksichtigung der Rundung. Ich verwende für Filter gerne auch 
statistische Rundung durch Rauschen, um Artefakte zu meiden. Um die 
Betrachtung des Rundens kommt man eh nicht herum, weil aufgrund der 
Rückkopplungen im IIR immer scheinbar wachsende Vektoren entstehen mit 
man Abschneiden muss. Ich frage mich nur immer noch, ob das bei dem CIC, 
der auch summiert oben korrekt gemacht ist.

[Gustl B. (-gb-)]

An welchen Stellen darf man denn bei einem IIR Filter abscheiden? Also 
wie im Eingangspost oben in dem Bildchen. Nur am Ausgang? Auch vor den 
einzelnen Registern? Wie ist das wenn man zwei dieser Biquads verbindet 
zu einer Kette, darf man dann zwischen diesen Beiden auch abschneiden?
Edit: Das müsste sogar egal sein weil bei so einer Kette keine 
Rückkopplung stattfindet.

FIR Filter sind mir irgendwie lieber ... oder eben real Zahlen.

[Marcel D. (diablokiller999)]

Jürgen S. schrieb:
> Ich frage mich nur immer noch, ob das bei dem CIC,
> der auch summiert oben korrekt gemacht ist.

Simulation und Test meiner Einheit im FPGA passen fast aufs Hz genau und 
der CIC sieht mir auch aus wie die normalen CIC Kurven. Denke, das 
sollte alles seine Richtigkeit haben.

Ich hab mir mein SOS-Filter nochmal mit anderen Koeffizienten gebastelt, 
dessen Polstellen nicht so nahe am Einheitskreis liegen und mit einem 
Butterworth-Verhalten klappt die Kaskadierung.

[J. S. (engineer)]

In dem anderen Bild weiter oben sah das etwas grob aus. Bei mir sehen 
Filter meistens so aus: http://96khz.org/doc/vasynthesisvhdl (Bild ganz 
unten).

[Marcel D. (diablokiller999)]

Mal noch eine Frage, da ich momentan Probleme mit dem DC Gain habe:
Bei einem Filter 4. Ordnung errechne ich den DC Gain für die 
SOS-Struktur ebenfalls mit dem Kehrwert der Division der Summe aller 
B-Koeffizienten durch die Summe aller A-Koeffizienten? Oder betrachte 
ich jeden Filter mit seinen Koeffizienten einzeln?

SOS-Filter 1
>b0=b2 = 1,2159e-004
>b1 = -1,3970e-004
>a1 = -1,8807e+000
>a2 = 8,8481e-001

SOS-Filter 2
>b0=b2 = 1,0000e+000
>b1 = -1,8295e+000
>a1 = -1,9450e+000
>a2 = 9,4941e-001

Wenn ich die Filter einzeln betrachte und den Gain-Faktor auf die 
einzelnen Filter anwende, kommt mein Ausgangssignal dem Eingangssignal 
zumindest nahe, aber so richtig zufrieden bin ich mit dem Ergebnis noch 
nicht.

[Martin O. (ossi-2)]

Gain(Hintereinanderschaltung)=Gain(Filter1) * Gain(Filter2)
gilt auch für DC, also Filter einzeln betrachten und DC-Verstärkungen 
multiplizieren.

Bei guter Skalierung und nicht zu kleinem Eingangssignal sollte sich der 
Gain=1 ziemlich genau ergeben.

Dein Filter1 hat ziemlich kleine bk Koeffizienten. Eventuell
verlierst Du da schon an Genauigkeit. Ich habe bei solchen
Filtern oft den DC-Gain sozusagen gleichmässig auf alle Stufen
verteilt. Das war gar nicht schlecht.

[Detlef _. (detlef_a)]

Ja, sehe ich so wie mein Vorredner. Die b des SOS1 zwei Zehnerpotenzen 
größer, die vom SOS2 2 Zehnerpotenzen kleiner. Ist aber nicht 
ungefährlich, weil SOS2 dann ein hundertmal größeres Eingangssignal 
sieht.

Cheers
Detlef

[Marcel D. (diablokiller999)]

Dann bin ich wohl zu blöd zum rechnen ^^
Ich habe mir jetzt ein neues Filter designed, die Koeffizienten sind:
SOS-Filter 1
>b0=b2 = 2,1754e-004
>b1 = 3,4375e-005
>a1 = -1,7550
>a2 = 7,7232e-001
> (b0+b1+b2)/(1+a1+a2) = 0,012811961189376
> Gain1 = 1/0,012811961189376 = 78,052062851174

SOS-Filter 2
>b0=b2 = 1,0000e+000
>b1 = -1,3305e+000
>a1 = -1,8827
>a2 = 9,0085e-001
> (b0+b1+b2)/(1+a1+a2) = 36,8870523415978
> Gain2 = 1/36,8870523415978 = 0,027109783420463

Wenn ich euch richtig verstanden habe, muss ich nun beide Gains 
multiplizieren
> Gges = Gain1 * Gain2 = 2,1159745194157

Alle B-Koeffizienten müssen nun mit diesem Gain multipliziert werden?
Das Filter ist laut PN-Plan stabil, ich habe die Bitbreite der 
Koeffizienten auf 30Bit erhöht. Leider bekomme ich bei meinem Sprung 
einen negativen Ausgangswert. Irgendwas kann da nicht stimmen, entweder 
ist meine Implementierung noch immer buggy (obwohl bei anderen Filtern 
die Filterkurve im FPGA stimmt) oder ich habe einen Denkfehler.

[Detlef _. (detlef_a)]

Nein, die a's spielen überhaupt keine Rolle, an denen kannst Du 
überhaupt nicht drehen. Es geht nur um die b's. Aus denen kannst Du 
einen beliebigen Koeffizienten rausziehen. Die b's von SOS1 sind grob 4 
Zehnerpotenzen kleiner als die von SOS2. Das haut auf die 
Rechengenauigkeit und die Festigkeit gegen Overflows. Also den einen 
Koeffezientensatz der b's 2 Zehnerpotenzen kleiner, den anderen 2 
Zehnerpotenzen grösser, hinten kommt dasgleiche raus.

Cheers
Detlef

[Marcel D. (diablokiller999)]

Detlef _. schrieb:
> Nein, die a's spielen überhaupt keine Rolle, an denen kannst Du
> überhaupt nicht drehen. Es geht nur um die b's.

An den a's habe ich auch nichts geändert ;)

Habe nun die B-Koeffizienten um 2 Bit verschoben, das zerhaut mir aber 
mein ganzes Filter. Mal die Guard Bits höher schrauben?

EDIT: Habe die Guardbits auf 8 hoch geschraubt, es verhält sich wieder 
wie ein Filter, aber mein Ausgangssignal ist noch immer 4x so hoch wie 
das Eingangssignal :-/

[Detlef _. (detlef_a)]

Allerhand schief:

- Wie kommst Du denn auf den b von 2.1754e-004? Das ist ungleich 
494253/(2^31) aus dem .png

- Skaliert QFORMAT 30 auf 2^30 oder aus 2^31 ? Vermutlich letzteres, 
weil sonst die a's keinen Sinn machen. Dann passen die b's aber 
überhaupt nicht.

- Die b's skalieren nur die Sprungantwort, die Form ändert sich nicht. 
Die simulierte Sprungantwort ist nicht richtig, zumindest nicht die, die 
Matlab richtig findet.

Overflow ist nach wie vor mein Tipp. Die Bitbreite hochziehen: 64+31 ~ 
100, zur Simulation mit intern 100Bit rechnen. Wenns denn geht kann man 
immer noch runterziehen. Der DC gain ist Dein kleinstes Problem, da kann 
man wenn alles fertig ist die b's skalieren. Das Ding klippt.

Cheers
Detlef

clear

fb1a=[2.1754e-004 3.4375e-005 2.1754e-004];
fa1a=[1 -1.7550 7.7232e-001];

fb2a=[1 -1.3305e+000 1 ];
fa2a=[1 -1.8827 9.0085e-001];

fb1=[494253/(2^31) 78100/(2^31) 494253/(2^31)];
fa1=[1 -1884416901/(2^30) 829272286/(2^30)];

fb2=[2272010340/(2^31) -3022909757/(2^31) 2272010340/(2^31)];
fa2=[1 -2021533732/(2^30) 967280322/(2^30)];


y=filter(fb1,fa1,ones(1,1000));
y=filter(fb2,fa2,y);
plot(y)
return

[Detlef _. (detlef_a)]

>>>>>>>>>>
EDIT: Habe die Guardbits auf 8 hoch geschraubt, es verhält sich wieder
wie ein Filter, aber mein Ausgangssignal ist noch immer 4x so hoch wie
das Eingangssignal :-/

Yo, weil Du die b's zweimal Faktor 2 falsch skalierst. Aber die 
Sprungantwort sieht ganz richtig aus.

Cheers
Detlef

[Marcel D. (diablokiller999)]

Detlef _. schrieb:
> Allerhand schief:
>
> - Wie kommst Du denn auf den b von 2.1754e-004? Das ist ungleich
> 494253/(2^31) aus dem .png
  (494253/((2^30)*Gges)) = 0,00021753994310336215479880413673356
  (78100/((2^30)*Gges)) = 0,000034374843564677572599026415780766
  (2272010340/((2^30)*Gges)) = 1
  (-3022909757/((2^30)*Gges)) = 1,3305000000897886654119436743203
Oder habe ich das mit dem Gesamt-Gain falsch verstanden?
> - Skaliert QFORMAT 30 auf 2^30 oder aus 2^31 ? Vermutlich letzteres,
> weil sonst die a's keinen Sinn machen.
Passt doch? => -1,755 * 2^30 = -1884416901,12

Ohne meine Gain-Korrektur lande ich ebenfalls bei 1

clear

fb1a=[2.1754e-004 3.4375e-005 2.1754e-004];
fa1a=[1 -1.7550 7.7232e-001];

fb2a=[1 -1.3305e+000 1 ];
fa2a=[1 -1.8827 9.0085e-001];

fb1=[233582/(2^30) 36910/(2^30) 233582/(2^30)];
fa1=[1 -1884416901/(2^30) 829272286/(2^30)];

fb2=[1073741824/(2^30) -1428613497/(2^30) 1073741824/(2^30)];
fa2=[1 -2021533732/(2^30) 967280322/(2^30)];


y=filter(fb1,fa1,ones(1,1000));
y=filter(fb2,fa2,y);
plot(y)

ya=filter(fb1a,fa1a,ones(1,1000));
ya=filter(fb2a,fa2a,ya);
plot(ya)

[Martin O. (ossi-2)]

Am besten erstmal unskaliertes und skaliertes Filter in C,Lua,Python 
oder
weiss Gott was simulieren und kontrollieren ob der DC-Gain dann
tatsächlich 1 ist.

[Detlef _. (detlef_a)]

Ja, die a's passen, die b's nicht. Gges ist die Gesamtverstärkung? Die 
verwurstest Du da nicht. Das ist ganz schlicht: Ist der Output Faktor 4 
zu hoch machst Du alle b's in SOS1 und SOS2 Faktor 2 kleiner. Der Output 
geht linear mit den b's mit.

Passt doch. Das zweite Filter hat auch noch eine Resonanzüberhöhung von 
ca. 3dB, entsprechend Faktor ~1.5. Dort muss Du auf jeden Fall noch 
prüfen ob er klippt, oder ohne auf die Resourcen zu achten die Guardbits 
noch eins hochziehen.

funzt.

Cheers
Detlef

[Marcel D. (diablokiller999)]

Martin O. schrieb:
> Am besten erstmal unskaliertes und skaliertes Filter in C,Lua,Python
> oder
> weiss Gott was simulieren und kontrollieren ob der DC-Gain dann
> tatsächlich 1 ist.
In Matlab geht es auf die 1 ohne Gain-Skalierung

Detlef _. schrieb:
>Ja, die a's passen, die b's nicht. Gges ist die Gesamtverstärkung? Die
>verwurstest Du da nicht.
Habe nun die Koeffizienten so genommen, wie ich sie aus MATLAB bekomme 
und nur mit 2^30 multipliziert. Gges besteht bei mir aus folgender 
Rechnung, habe ich aber in der jetzigen Simulation außen vor gelassen um 
jeden Rechenfehler auszuschließen:

SOS-Filter 1
>b0=b2 = 2,1754e-004
>b1 = 3,4375e-005
>a1 = -1,7550
>a2 = 7,7232e-001
> (b0+b1+b2)/(1+a1+a2) = 0,012811961189376
> Gain1 = 1/0,012811961189376 = 78,052062851174

SOS-Filter 2
>b0=b2 = 1,0000e+000
>b1 = -1,3305e+000
>a1 = -1,8827
>a2 = 9,0085e-001
> (b0+b1+b2)/(1+a1+a2) = 36,8870523415978
> Gain2 = 1/36,8870523415978 = 0,027109783420463
> Gges = Gain1 * Gain2 = 2,1159745194157

=> Jetzt alle b's mit dem Gges multiplizieren?

Nun bekomme ich aber mit den blanken Koeffizienten einen negativen 
Ausgang.

EDIT: In der Matlab Simulation waren die B's nicht 1 sondern 1073741824, 
falsche Beschriftung aber Simulation lief mit den richtigen Werten :D

[Martin O. (ossi-2)]

>>In Matlab geht es auf die 1 ohne Gain-Skalierung

Wenn Du in Matlab und VHDL die gleichen Koeffizienten mit den gleichen 
Skalierungen verwendest sollte eigentlich bei beiden auch der gleiche 
Gain rauskommen.

[Martin O. (ossi-2)]

@ Marcel D.
Ich würde an Deiner Stelle in Deine VHDL Simulation mal numerische 
Testausgaben einbauen, damit Du genau nachvollziehen kannst, was 
passiert. Und dann wird solange nach Fehlern gesucht, bis Matlab und 
VHDL aufs letzte Bit übereinstimmen. Alles andere ist Fischen und 
Vermuten im Trüben.

[Detlef _. (detlef_a)]

Marcel, machs nicht so umständlich. Es war Faktor 4 zu gross, dann machs 
zweimal Faktor 2 kleiner und fertig ist die Laube.

Bei Deinen neuen Versuchen hast Du vllt. wieder nen klipping bei den b's 
eingebaut. Wenn die signed sind und der Wert zu gross wird kriegst Du 
einen negativen Ausgang.

Cheers
Detlef

[Marcel D. (diablokiller999)]

Detlef _. schrieb:
> Marcel, machs nicht so umständlich. Es war Faktor 4 zu gross, dann machs
> zweimal Faktor 2 kleiner und fertig ist die Laube.
>
So, es ist vollbracht, habe die Koeffizienten um jeweils zwei 
Zehnerpotenzen nach links/rechts verschoben, meine neuen Koeffizienten 
sehen nun wie folgt aus:
SOS-Filter 1
>b0=b2 = 0,021754
>b1 = 0,0034375
>a1 = -1,7550
>a2 = 7,7232e-001

SOS-Filter 2
>b0=b2 = 0,01
>b1 = -0,013305
>a1 = -1,8827
>a2 = 9,0085e-001

Auf 2^30 skaliert, bekomme ich ein annähernd ähnliches Ausgangs- wie 
Eingangssignal (siehe Simulation). Jedoch fangen die Werte ab einem 
Eingangswert über 0FFFFFFFFFFFFFFF wieder an, auseinander zu driften 
(Ausgang < Eingang). Dabei spielt die Anzahl der Guard Bits keine Rolle, 
selbst auf 16 ändert sich nichts am Verhalten.

[Mitleser (Gast)]

und ein weiteres Aufdrehen der Bits willst Du nicht?
was tut sich da?

[Detlef _. (detlef_a)]

Die Eingangs/Ausgangs Bitbreiten aufdrehen, sehe ich auch so, das ist 
sicher Eingangs/Ausgangsclipping. Wenn Du mit 60Bit Werten reingehst 
hast Du noch 3Bit Luft, das ist nicht viel. Die DC gains des SOS1 und 
SOS2 sind Faktor ~2 unterschiedlich, da ist eines der Bits schon weg.

Cheers
Detlef

Und ach ja: Bei einem des SOS hatte ich ja den Faktor 1.5 als 
Resonanzüberhöhung berechnet, da ist das zweite Bit perdu. Das ist alles 
sehr auf Kante genäht, da darf man nicht viele Fehler machen. Jetzt mal 
bißchen Großzügigkeit zeigen ;) ( ausserdem brauchst du keine 64Bit)

[Edi M. (Gast)]

Detlef _. schrieb:
> Wenn Du mit 60Bit Werten reingehst

dann hat er aber was falsch gemacht. Kein Mensch brauch irgendwo 60 
Bit-Auflösung.

[Marcel D. (diablokiller999)]

Edi M. schrieb:
> Detlef _. schrieb:
>> Wenn Du mit 60Bit Werten reingehst
>
> dann hat er aber was falsch gemacht. Kein Mensch brauch irgendwo 60
> Bit-Auflösung.

Ich hab das Teil nicht designed, ich soll es nur bauen ^^
Soweit ich es sehe, sollte meine FPGA Implementierung auch richtig sein. 
Zumindest verhält sich das Filter als einzelnes Filter 2. Ordnung so, 
wie es soll. Eventuell ist meine Filterkurve zu sehr auf Kante genäht 
und bei Werten über 60Bit wird es halt instabil?

EDIT:
Ich habe mir mal einen Butterworth Filter 4. Ordnung designed und die 
Koeffizienten übertragen. Das Filter scheint wunderbar zu funktionieren, 
zumindest wenn ich Sprünge drauf gebe ist der Ausgangswert nahezu 
identisch mit dem Eingangswert. Liegt also wirklich eher an den 
Filterkoeffizienten wie es scheint.

[Detlef _. (detlef_a)]

>>Liegt also wirklich eher an den Filterkoeffizienten wie es scheint.

Ja, das liegt an den Filtereigenschaften. Die Sprungantwort des Filters 
schwingt über, also hast Du auf jeden Fall clipping wenn Du mit deinem 
Maximalsignal reingehst. Wieviel kleiner als das Maximalsignal das 
Eingangssignal sein darf, ohne dass es zu clipping kommt hängt von der 
Höhe des Überschwingens ab.

>>>
Soweit ich es sehe, sollte meine FPGA Implementierung auch richtig sein.

Ja, das Filter ist ok, sehe ich auch so.

>>>
Eventuell ist meine Filterkurve zu sehr auf Kante genäht und bei Werten 
über 60Bit wird es halt instabil?

Nein, Stabilität ist eine Eigenschaft des Filters und hat mit dem 
Eingangssignal nix zum tun. Das Ding clippt, weil die Ausgangsbitbreiten 
bei hoch ausgesteuertem Eingang zu klein sind, der Ausgang braucht 
'Platz zum Überschwingen'.

Mit genau dieser Argumentation erzählst Du Deinem 'Auftraggeber', dass 
er nicht mit 64 voll ausgesteuerten Bit reingehen kann. Wenn er das 
will, muss die Eingangs/Ausgangsbitbreite grösser werden. Wenn Du Lust 
und die Möglichkeit zum bashen hast, kannst Du deinem Auftraggeber 
gleich auch noch erzählen, dass die vorgegebene Spec. im wesentlichen 
sein eigenes technisches Unverständnis dokumentiert.

Schönes W'ende.
Cheers
Detlef

[Martin O. (ossi-2)]

@ Marcel D.
Könntest Du Deinen aktuellen VHDL Code nochmal hier reinstellen
(inkl. aktueller Koeffizienten)? Ich würde da gerne nochmal was 
ausprobieren.

[Michael W. (Gast)]

Martin O. schrieb:
> @ Marcel D.
> Könntest Du Deinen aktuellen VHDL Code nochmal hier reinstellen
> (inkl. aktueller Koeffizienten)? Ich würde da gerne nochmal was
> ausprobieren.

Hausaufgaben saugen?

[Martin O. (ossi-2)]

Ich habe jetzt JAVA-Routinen geschrieben, mit denen man die
Arithmetik solcher Filter simulieren kann. Dabei wird
auch automatisch überwacht, ob irgendwo ein Überlauf stattfindet.
Gleichzeitig wird für jeden Wert der maximal aufgetretene Wert 
ermittelt.

@ Marcel D
Damit könnte ich Dein Filter mal analysieren.

@ elektrowagi78
Bevor Du solche Unterstellungen loslässt, lies Dir lieber mal
den Thread durch, und entscheide dann, ob ich eher Geber oder Nehmer
in diesem Thread bin.

[Edi M. (Gast)]

Marcel D. schrieb:
> Ich hab das Teil nicht designed, ich soll es nur bauen ^^

Dann frage Ich doch mal anders herum: Wer hat es womit DESIGNED?
Und mit welcher Methode und Richtlinie wurde es (automatisch) übersetzt?

Martin O. schrieb:
> Gleichzeitig wird für jeden Wert der maximal aufgetretene Wert
> ermittelt.
Ich nehme ab, für jede Veriable wird der maximale Wert ermittelt. Wie 
testest Du das Szenatio? Wie ermittelst Du den testvector für den worst 
case, den Du auf den Filter loslassen musst, damit Du den worst case 
auch erzeugst?

[Marcel D. (diablokiller999)]

Martin O. schrieb:
> @ Marcel D.
> Könntest Du Deinen aktuellen VHDL Code nochmal hier reinstellen
> (inkl. aktueller Koeffizienten)? Ich würde da gerne nochmal was
> ausprobieren.

Koeffizienten für das erste Filter der Kaskade:
>A1 => -1465764964,
>A2 => 512819095,
>B0 => 30199151,
>B1 => 60397654,
>B2 => 30199151,
>GUARDBITS => 4,
>INPUT_WIDTH => 64,
>OUTPUT_WIDTH => 64,
>QFORMAT => 30

Koeffizienten für das zweite Filter der Kaskade:
>A1 => -1730549698,
>A2 => 799411525,
>B0 => 35648239,
>B1 => 71307173,
>B2 => 35648239,
>GUARDBITS => 4,
>INPUT_WIDTH => 64,
>OUTPUT_WIDTH => 64,
>QFORMAT => 30

@Edi:
Die Designvorgaben waren, möglichst steiles Filterverhalten bei 
möglichst wenig Platz auf dem FPGA und 64Bit Datenbreite. Habe mich nun 
in Matlab eingearbeitet und ein paar Filter durchprobiert, bis ich auf 
Koeffizienten gestoßen bin, die meiner vagen Aufgabenbeschreibung 
entsprechen. Mit den aktuellen Koeffizienten habe ich laut Simulation 
auch eine Verstärkung von 1, leider weicht das ganze in der Realität nun 
extrem davon ab. Ersetze ich das Filter durch ein FIR-Design, stimmt die 
Amplitude wieder.

[Martin O. (ossi-2)]

@ Marcel D
Kannst Du mir auch noch die Fliesspunktwerte der Filterkoeffizienten 
geben?

[Marcel D. (diablokiller999)]

Martin O. schrieb:
> @ Marcel D
> Kannst Du mir auch noch die Fliesspunktwerte der Filterkoeffizienten
> geben?
Filter1:
A1 = -1,37E+00
A2 = 4,78E-01
B0 = 9,33E-04
B1 = 1,87E-03
B2 = 9,33E-04

Filter2:
A1 = -1,61E+00
A2 = 7,45E-01
B0 = 1
B1 = 2,00E+00
B2 = 1

Die Koeffizienten für die FPGA Implementierung sind korrigiert, um einen 
DCGain von 1 zu erhalten. Also nicht wundern, wenn du mit 2^30 
multiplizierst auf andere Werte kommst ;)

[Martin O. (ossi-2)]

Ok, ich weiss ich nerve.....
Kannst Du mir die Fliespunkt-Koeffizienten mal mit mindestens 6 Stellen
Genauigkeit geben. Genauso die Normalöieesierungskoeffzzienten, die
Du nimmst, um im FPGA gain=1 zu erreichen.

[Marcel D. (diablokiller999)]

Martin O. schrieb:

> Kannst Du mir die Fliespunkt-Koeffizienten mal mit mindestens 6 Stellen
> Genauigkeit geben.
Gibt mir MATLAB genau so aus, genauer kann ich sie dir nicht liefern
> Genauso die Normalöieesierungskoeffzzienten, die
> Du nimmst, um im FPGA gain=1 zu erreichen.

Ich addiere alle B's pro Stufe und dividiere sie durch die Summe aller 
A's.
Komme für die erste Stufe auf 1/0,033187733 und für die Zweite auf 
1/30,12047286. Damit multipliziere ich die einzelnen Koeffizienten des 
jeweiligen Filters und erhalte in der Simulation (fast genau) die 
gleiche Amplitude wie am Eingang.

[Martin O. (ossi-2)]

Die Werte von Filter1 kann ich jetzt nachvollziehen,
aber bei Filter2 geht was schief:

A1=-1.61
A1*2^30=-1.61*1073741824=-1728724337
Du gibst für VHDL aber >A1 => -1730549698,

[Marcel D. (diablokiller999)]

Martin O. schrieb:

> A1=-1.61
> A1*2^30=-1.61*1073741824=-1728724337
> Du gibst für VHDL aber >A1 => -1730549698,

Hab doch genauere Werte bekommen xD
Filter1:
A1 = -1,3651
A2 = 0,4776
B0 = 0,00093341
B1 = 0,0018668
B0 = 0,00093341
Gain = 30,1316148938564

Filter2:
A1 = -1,6117
A2 = 0,74451
B0 = 1
B1 = 2,0003
B0 = 1
Gain = 0,03320000999925

[Martin O. (ossi-2)]

So, eine erste Version von Routinen zur Simulation von FPGA Arithmetik/Filtern in Java ist fertig.

Werte werden in "Variablen" vom Typ "DSPvalue" gespeichert. In "DSPvalue"-Werten steht nicht
nur der eigentliche Variablenwert (als Java "BigInteger"), sondern auch Skalierungsfaktor, 
Anzahl der Bits zur Speicherung  und einiges mehr. Man kann für jede "DSPvalue" Variable
ein "logging" einschalten, dann wird jede Zuweisung protokolliert.
Eine Variable mit "varWidth" Bits, skaliert mit "scale", Anfangswert=0 und Namen "IIRoutput1"
erzeugt man mittels:

DSPvalue IIRoutput1=new DSPvalue( 0 ,varWidth ,scale1,"IIRoutput1" ) ;

Zum Rechnen gibt es ein "DSP"-Objekt, damit kann man die Standardoperationen auf Variablen
vom Typ "DSPvalue" ausführen. Dabei wird überwacht, ob ein Überlauf stattfindet. Weiter
kann man eine Protokollierung einschalten, dann werden alle Operationen ausgegeben.

Der Code, der das IIR Filter von Marcel D. simuliert, sieht wie folgt aus (vergleiche VHDL Code von Marcel):  
  
  
void step(long IIRinput, DSPvalue IIRout){
  nX.valueSet(IIRinput)  ;
  
//VHDL:                        nA1 <= nZ1 * cA1;   
  DSP.mul(nZ1,cA1,nA1) ;
//VHDL:                        nA2 <= nZ2 * cA2; 
  DSP.mul(nZ2,cA2,nA2) ;
//VHDL:                        nB1 <= nZ1 * cB1;  
  DSP.mul(nZ1,cB1,nB1) ;
//VHDL:                        nB2 <= nZ2 * cB2;  
  DSP.mul(nZ2,cB2,nB2) ;
//VHDL:                        nX  <= shift_left(signed(nX),QFORMAT);   
  DSP.shl(nX,scaleShift,nXhelp) ;
//VHDL:                        nSUMA1A2 <= nA1(nA1'left-1 downto 0) + nA2(nA2'left-1 downto 0);   
  DSP.add(nA1,nA2,nSUMA1A2) ;
  DSP.add(nB1,nB2,nSUMB1B2) ;
  DSP.sub(nXhelp,nSUMA1A2,nSUMx) ;

  DSP.shr(nSUMx,scaleShift,nSUMxHelp) ;
  
  DSP.copy(nSUMxHelp,nZ1help) ;
  DSP.copy(nZ1,nZ2help ) ;
  DSP.mul(nSUMxHelp,cB0,nB0) ;
        
  DSP.copy(nZ1help,nZ1 ) ;
  DSP.copy(nZ2help,nZ2 ) ;
    
  DSP.add(nB0,nSUMB1B2,nSUMxB0) ;
  
  DSP.shr(nSUMxB0,scaleShift,IIRout) ;
  }
    
Der Code zur Simulation zweier IIR Filters hintereinander sieht dann wie folgt aus
(Eingangssprung von  0 auf 10000, Filter2 schwingt etwas über):
  
nSamples=100 ;  
for(int k=0 ; k<nSamples ; k++){
  if(k>=10) { IIRinput=10000 ; } else { IIRinput=0 ; }
  IIRfil1.step(IIRinput,IIRoutput1) ;
  IIRfil2.step(IIRoutput1.asBigInt.longValue(),IIRoutput2) ;
  }  
  
Die Ausgabe ergibt sich dann als :
  
  Fil1.nX :=       0 IIRoutput1 :=       0 IIRoutput2 :=       0
  Fil1.nX :=       0 IIRoutput1 :=       0 IIRoutput2 :=       0
  Fil1.nX :=   10000 IIRoutput1 :=     281 IIRoutput2 :=       9
  Fil1.nX :=   10000 IIRoutput1 :=    1227 IIRoutput2 :=      74
  Fil1.nX :=   10000 IIRoutput1 :=    2666 IIRoutput2 :=     292
  Fil1.nX :=   10000 IIRoutput1 :=    4178 IIRoutput2 :=     771
  Fil1.nX :=   10000 IIRoutput1 :=    5555 IIRoutput2 :=    1576
  Fil1.nX :=   10000 IIRoutput1 :=    6713 IIRoutput2 :=    2697
  Fil1.nX :=   10000 IIRoutput1 :=    7635 IIRoutput2 :=    4056
  Fil1.nX :=   10000 IIRoutput1 :=    8342 IIRoutput2 :=    5536
  Fil1.nX :=   10000 IIRoutput1 :=    8866 IIRoutput2 :=    7005
  Fil1.nX :=   10000 IIRoutput1 :=    9243 IIRoutput2 :=    8340
  Fil1.nX :=   10000 IIRoutput1 :=    9509 IIRoutput2 :=    9451
  Fil1.nX :=   10000 IIRoutput1 :=    9691 IIRoutput2 :=   10282
  Fil1.nX :=   10000 IIRoutput1 :=    9812 IIRoutput2 :=   10821
  Fil1.nX :=   10000 IIRoutput1 :=    9891 IIRoutput2 :=   11086
  Fil1.nX :=   10000 IIRoutput1 :=    9941 IIRoutput2 :=   11124
  Fil1.nX :=   10000 IIRoutput1 :=    9972 IIRoutput2 :=   10994
  Fil1.nX :=   10000 IIRoutput1 :=    9989 IIRoutput2 :=   10761
  ......
  Fil1.nX :=   10000 IIRoutput1 :=    9999 IIRoutput2 :=    9998
  Fil1.nX :=   10000 IIRoutput1 :=    9999 IIRoutput2 :=    9998
  Fil1.nX :=   10000 IIRoutput1 :=    9999 IIRoutput2 :=    9998

Name= Fil1.nA1            absMax/scale=   121411.99392  absMax=  130365135809220.00 [  46.9 Bits]  of [48] Bits
Name= Fil1.nA2            absMax/scale=    42477.74399  absMax=   45610130309300.00 [  45.4 Bits]  of [48] Bits
Name= Fil1.nB0            absMax/scale=     2501.45085  absMax=    2685912401000.00 [  41.3 Bits]  of [48] Bits
Name= Fil1.nB1            absMax/scale=     5002.84820  absMax=    5371767346760.00 [  42.3 Bits]  of [48] Bits
Name= Fil1.nB2            absMax/scale=     2501.45085  absMax=    2685912401000.00 [  41.3 Bits]  of [48] Bits
Name= Fil1.nZ1            absMax/scale=    88940.00000  absMax=            88940.00 [  16.4 Bits]  of [48] Bits
Name= Fil1.nZ2            absMax/scale=    88940.00000  absMax=            88940.00 [  16.4 Bits]  of [48] Bits
Name= Fil1.nSUMA1A2       absMax/scale=    78940.79353  absMax=   84762031636924.00 [  46.3 Bits]  of [48] Bits
Name= Fil1.nSUMB1B2       absMax/scale=     7504.24280  absMax=    8057619350106.00 [  42.9 Bits]  of [48] Bits
Name= Fil1.nSUMx          absMax/scale=    88940.79353  absMax=   95499449876924.00 [  46.4 Bits]  of [48] Bits
Name= Fil1.nSUMxHelp      absMax/scale=    88940.00000  absMax=            88940.00 [  16.4 Bits]  of [48] Bits
Name= Fil1.nSUMxB0        absMax/scale=    10005.46865  absMax=   10743290157906.00 [  43.3 Bits]  of [48] Bits
Name= Fil1.nZ1help        absMax/scale=    88940.00000  absMax=            88940.00 [  16.4 Bits]  of [48] Bits
Name= Fil1.nZ2help        absMax/scale=    88940.00000  absMax=            88940.00 [  16.4 Bits]  of [48] Bits
Name= Fil1.nX             absMax/scale=    10000.00000  absMax=            10000.00 [  13.3 Bits]  of [48] Bits
Name= Fil1.nXhelp         absMax/scale=    10000.00000  absMax=   10737418240000.00 [  43.3 Bits]  of [48] Bits
Name= Fil1-cB0            absMax/scale=        0.02813  absMax=         30199150.00 [  24.8 Bits]  of [32] Bits
Name= Fil1-cB1            absMax/scale=        0.05625  absMax=         60397654.00 [  25.8 Bits]  of [32] Bits
Name= Fil1-cB2            absMax/scale=        0.02813  absMax=         30199150.00 [  24.8 Bits]  of [32] Bits
Name= Fil1-cA1            absMax/scale=       -1.36510  absMax=      -1465764963.00 [  30.4 Bits]  of [32] Bits
Name= Fil1-cA2            absMax/scale=        0.47760  absMax=        512819095.00 [  28.9 Bits]  of [32] Bits
Name= IIRoutput1          absMax/scale=    10005.00000  absMax=            10005.00 [  13.3 Bits]  of [48] Bits
Name= Fil2.nA1            absMax/scale=   135510.12424  absMax=  145502887974063.00 [  47.0 Bits]  of [48] Bits
Name= Fil2.nA2            absMax/scale=    62597.65626  absMax=   67213721610475.00 [  45.9 Bits]  of [48] Bits
Name= Fil2.nB0            absMax/scale=     2791.42362  absMax=    2997268286881.00 [  41.4 Bits]  of [48] Bits
Name= Fil2.nB1            absMax/scale=     5583.68470  absMax=    5995435798667.00 [  42.4 Bits]  of [48] Bits
Name= Fil2.nB2            absMax/scale=     2791.42362  absMax=    2997268286881.00 [  41.4 Bits]  of [48] Bits
Name= Fil2.nZ1            absMax/scale=    84079.00000  absMax=            84079.00 [  16.4 Bits]  of [48] Bits
Name= Fil2.nZ2            absMax/scale=    84079.00000  absMax=            84079.00 [  16.4 Bits]  of [48] Bits
Name= Fil2.nSUMA1A2       absMax/scale=    74188.17233  absMax=   79658943478564.00 [  46.2 Bits]  of [48] Bits
Name= Fil2.nSUMB1B2       absMax/scale=     8369.56392  absMax=    8986750829635.00 [  43.0 Bits]  of [48] Bits
Name= Fil2.nSUMx          absMax/scale=    84079.17233  absMax=   90279323859748.00 [  46.4 Bits]  of [48] Bits
Name= Fil2.nSUMxHelp      absMax/scale=    84079.00000  absMax=            84079.00 [  16.4 Bits]  of [48] Bits
Name= Fil2.nSUMxB0        absMax/scale=    11124.40113  absMax=   11944734757138.00 [  43.4 Bits]  of [48] Bits
Name= Fil2.nZ1help        absMax/scale=    84079.00000  absMax=            84079.00 [  16.4 Bits]  of [48] Bits
Name= Fil2.nZ2help        absMax/scale=    84079.00000  absMax=            84079.00 [  16.4 Bits]  of [48] Bits
Name= Fil2.nX             absMax/scale=    10005.00000  absMax=            10005.00 [  13.3 Bits]  of [48] Bits
Name= Fil2.nXhelp         absMax/scale=    10005.00000  absMax=   10742786949120.00 [  43.3 Bits]  of [48] Bits
Name= Fil2-cB0            absMax/scale=        0.03320  absMax=         35648239.00 [  25.1 Bits]  of [32] Bits
Name= Fil2-cB1            absMax/scale=        0.06641  absMax=         71307173.00 [  26.1 Bits]  of [32] Bits
Name= Fil2-cB2            absMax/scale=        0.03320  absMax=         35648239.00 [  25.1 Bits]  of [32] Bits
Name= Fil2-cA1            absMax/scale=       -1.61170  absMax=      -1730549697.00 [  30.7 Bits]  of [32] Bits
Name= Fil2-cA2            absMax/scale=        0.74451  absMax=        799411525.00 [  29.6 Bits]  of [32] Bits
Name= IIRoutput2          absMax/scale=    11124.00000  absMax=            11124.00 [  13.4 Bits]  of [48] Bits

Man erkennt, das in beiden Filtern der nA1 Wert am grössten ist, un d deshalb am ehesten einen Überlauf
verursachen wird. Hier reichen die vorgesehenen Bits (48) gertade aus, um die nötigen Bits (46.9) zu speichern.
Viel Luft ist da nicht mehr, was man bei einem größeren Eingangssignal auch merkt: 


Bei einem Sprung von 0 auf 20000 am Eingang ergibt sich ein Überlauf bei Wert Fil2.nA1:
    
Error: MUL overflow 2
Fil2.nA1 =mul(Fil2.nZ1=164017 * Fil2-cA1=-1730549697 ) 
Name= Fil2.nA1  nBits= 48  Value=-269957099983515 [  47.9 Bits] value/scale=-251417.141392375340  
asDouble=-269966856015344.06000  range is -140737488355328 .. 140737488355327   [  47.0 Bits] 
  
Statt eines Eingangssprungs kann man natürlich auch einen Frequenzsweep oder andere Eingangssignale 
simulieren. Wie man sieht kann man mit diesen Routinen gut untersuchen, was eine Filterimplementation so macht.

[Martin O. (ossi-2)]

Hab vergessen: Diese Tabelle gibt Auskunft darüber, wie gross die
entsprechende Variable maximal wurde, wievuiel Bits sie dann braucht, 
und wieviele Bits reserviert sind

Name= Fil1.nA1  absMax/scale=   121411.99392  absMax=  130365135809220.00 [  46.9 Bits]  of [48] Bits
Name= Fil1.nA2  absMax/scale=    42477.74399  absMax=   45610130309300.00 [  45.4 Bits]  of [48] Bits

Die Variable Fil1.nA1 ist kurz davor (46.9 Bits von 48 Bits) einen
Überlauf zu verursachen.

[W.S. (Gast)]

Ich bin erschlagen!
Wie benutzt man das?

[Martin O. (ossi-2)]

Hier mal ein "vollständiges" Beispiel. Zuerst der JAVA Code:

package allClasses;

public class Demo1v01 {
/*  A1=-1.0*scale
  A2= 0.5*scale
  B0= 1.0*scale
  yNew=(B0*x0-A1*y1-A2*y2)/scale
  y2=y1 ;
    y1=yNew
*/
  static void demo(){            
    int scaleShift=16 ;
    long scale=1L<<scaleShift ;
    int scale1=1 ;
    int coeffWidth=17 ;
    int varWidth=16 ;
    int productWidth=32 ;
    
    DSPvalue B0=new DSPvalue( (long)( 1.0*scale) ,coeffWidth ,scale,"B0" ) ;
    DSPvalue A1=new DSPvalue( (long)(-1.0*scale) ,coeffWidth ,scale,"A1" ) ;
    DSPvalue A2=new DSPvalue( (long)( 0.5*scale) ,coeffWidth ,scale,"A2" ) ;
    
    DSPvalue x0=new DSPvalue( 0 ,varWidth ,scale1,"x0" ) ;
    DSPvalue y1=new DSPvalue( 0 ,varWidth ,scale1,"y1" ) ;
    DSPvalue y2=new DSPvalue( 0 ,varWidth ,scale1,"y2" ) ;
    DSPvalue yNew=new DSPvalue( 0 ,varWidth ,scale1,"yNew" ) ;
  
    DSPvalue pB0=new DSPvalue( 0 ,productWidth ,scale,"pB0" ) ;
    DSPvalue pA1=new DSPvalue( 0 ,productWidth ,scale,"pA1" ) ;
    DSPvalue pA2=new DSPvalue( 0 ,productWidth ,scale,"pA2" ) ;
    DSPvalue sum=new DSPvalue( 0 ,productWidth ,scale,"sum" ) ;
  
    yNew.valueProt=true ;
    x0.valueProt=true ;
    pA1.valueProt=true ;
    pA2.valueProt=true ;
  
    int nSamples=10 ;
    int  IIRinput ;
    for(int k=0 ; k<nSamples ; k++){
      if(k==2) { IIRinput=10000 ; } else { IIRinput=0 ; }
      x0.valueSet(IIRinput) ;
      DSP.mul(x0, B0, pB0) ;         // pB0 = x0 * B0
      DSP.mul(y1, A1, pA1) ;         // pA1 = y1 * A1
      DSP.mul(y2, A2, pA2) ;         // pA2 = y2 * A2
      DSP.copy(pB0,sum) ;            // sum = pB0
      DSP.sub(sum,pA1,sum) ;         // sum=pB0-pA1
      DSP.sub(sum,pA2,sum) ;         // sum=pB0-pA1-pA2
      DSP.shr(sum,scaleShift,yNew) ; // yNew=sum>>scaleShift
      DSP.copy(y1,y2) ;              // y2 = y1
      DSP.copy(yNew,y1) ;            // y1=yNew
      System.out.println() ;
      }
    DSPvalue.maxListShow() ;  
    }
  }

Das Programm erzeugt die folgende Ausgabe:

       x0 :=         0      pA1 :=         0      pA2 :=         0     yNew :=         0
       x0 :=         0      pA1 :=         0      pA2 :=         0     yNew :=         0
       x0 :=     10000      pA1 :=         0      pA2 :=         0     yNew :=     10000
       x0 :=         0      pA1 :=-655360000      pA2 :=         0     yNew :=     10000
       x0 :=         0      pA1 :=-655360000      pA2 := 327680000     yNew :=      5000
       x0 :=         0      pA1 :=-327680000      pA2 := 327680000     yNew :=         0
       x0 :=         0      pA1 :=         0      pA2 := 163840000     yNew :=     -2500
       x0 :=         0      pA1 := 163840000      pA2 :=         0     yNew :=     -2500
       x0 :=         0      pA1 := 163840000      pA2 := -81920000     yNew :=     -1250
       x0 :=         0      pA1 :=  81920000      pA2 := -81920000     yNew :=         0

Das sind die Daten der Filtersimulation. Danach wird
die folgende Tabelle ausgegeben, die für jede Variable
den Maximalwert ausgibt

Name= B0        absMax/scale=        1.00000  absMax=       65536.00 [  16.0 Bits]  of [17] Bits
Name= A1        absMax/scale=       -1.00000  absMax=      -65536.00 [  16.0 Bits]  of [17] Bits
Name= A2        absMax/scale=        0.50000  absMax=       32768.00 [  15.0 Bits]  of [17] Bits
Name= x0        absMax/scale=    10000.00000  absMax=       10000.00 [  13.3 Bits]  of [16] Bits
Name= y1        absMax/scale=    10000.00000  absMax=       10000.00 [  13.3 Bits]  of [16] Bits
Name= y2        absMax/scale=    10000.00000  absMax=       10000.00 [  13.3 Bits]  of [16] Bits
Name= yNew      absMax/scale=    10000.00000  absMax=       10000.00 [  13.3 Bits]  of [16] Bits
Name= pB0       absMax/scale=    10000.00000  absMax=   655360000.00 [  29.3 Bits]  of [32] Bits
Name= pA1       absMax/scale=    10000.00000  absMax=   655360000.00 [  29.3 Bits]  of [32] Bits
Name= pA2       absMax/scale=     5000.00000  absMax=   327680000.00 [  28.3 Bits]  of [32] Bits
Name= sum       absMax/scale=    10000.00000  absMax=   655360000.00 [  29.3 Bits]  of [32] Bits
Stop :JavaTestBed1v01