Forum: Digitale Signalverarbeitung / DSP / Machine Learning iir Filter als Sinusgenerator


von Klausi (Gast)


Lesenswert?

Hallo,
habe mir nach einer Vorlage aus'm Netz einen iir-Filter 2. Ordnung 
gebaut und der schwingt auch schön:

Formel:
y[n]= A*y[n-1] + B*y[n-2] + x[n]
y[n]= 1,9757*y[n-1] - y[n-2] + x[n]

bestimmt mit:

y[0] = sin((2pi/40)*0) = 0
y[1] = sin((2pi/40)*1) = 0,1564
y[2] = sin((2pi/40)*2) = 0,3090
y[2] = sin((2pi/40)*3) = 0,0454

eingesetzt in:
y[2] = A*y[1] + B*y[0]
y[3] = A*y[2] + B*y[1]

findet man:
B = -1
A = 1,9757

Formel:
y[n]= A*y[n-1] + B*y[n-2] + x[n]
y[n]= 1,9757*y[n-1] - y[n-2] + x[n]

alles in "int 32bit" gepackt und in eine Funktion geschrieben kann ich 
die Schwingung hören.

Probleme habe ich mit der Bestimmung der Frequenz:
Kann man die Aufteilung einer Schwingung(2pi) in 40 Teile als 
"Schrittweite" interpretieren? -> Nach 40 Aufrufen der iir() Funktion 
erhalte ich eine komplette Sinus-Schwingung.

Ich rufe die Funktion mit ca 48 kHz auf, somit erhalte ich
T(Sinusschwingung) = 40 * 1/48 kHz = 0,8333 ms entspricht ca. 1200 Hz.
Stimmt das soweit? Wie kann ich die Schwingung "feiner" machen bzw. die 
Frequenz verändern.

Ich habe schon den Koeffizienten A für 2pi/1meg ausgegerchnet aber dann 
schwingt der Filter nicht mehr so schön rund.

Frequenz verändern durch varieren des Aufrufintervalls (z.B. vergrößern) 
der Funktion macht die Schwingung auch immer "gröber".

Kann mir jemand sagen, ob man die Frequenz der Schwingung des Filters 
dem iir-Filter praktisch implizit beibringen kann?

Viele Grüße

von J. S. (engineer) Benutzerseite


Lesenswert?

Das Filter hat eine Grenzfrequenz, die Du berechnen und auf das Filter 
anpassen kannst. Die anregende Frequenz muss dazu passen.

"nicht so schön rund" deutet darauf hin, dass die Anregung zu schnell 
ist, bzw die Grenzfrequenz nicht nachgeregelt ist.

Am Besten spielst Du Dir das Filterverhalten in Excel oder SciCos / 
Matlab mal durch.

1. Spalte: alle 10 Werte alternierende 100, -100
2. Spalte: Wert = 0,1 * erste Spalte + 0,9 * alter Wert obendrüber
3. Spalte: Wert = 0,1 * zweite Spalte + 0,9 * alter Wert obendrüber

Schon ab der 2. Stufe macht das aus einem Rechteck einen reltaiv guten 
Sinux mit -6dB Amplitude.

Für geringere Frequenzen ist die Verteilung nue zu alt dann z.B. 0,05 / 
0,95

Aus der Anstiegszeit des Filters kann man die Grenzfrequenz bestimmen. 
Der Faktor war glaube ich e, jedenfalls knapp Faktor 3.  f = 2,8 * 1/t

Ich hänge so einen Filter gerne hinter meine DDSen, wenn ich mit kleinen 
Tabellen arbeite. Man muss nur Phase und Amplitude korrigieren.

von Christoph db1uq K. (christoph_kessler)


Angehängte Dateien:

Lesenswert?

Das Buch ISBN 9780792395591 habe ich mir gerade gebraucht gekauft, darin 
ist auch ein IIR-Sinus/Cosinus-Generator beschrieben. Anscheinend sollte 
man spätestens nach ein paar Schwingungen wieder auf die 
Startbedingungen  resetten.

von Detlef _. (detlef_a)


Lesenswert?

Der A Koeffizient bestimmt Dir die Frequenz, B ist immer -1 .

Rechnung geht so:
1.9757 = 2*cos(Schrittweite)
dann ist Schrittweite = 0.1560 rad
In 2*pi/Schrittweite = 40.265 Schritten bist Du also einmal 'rum' um den 
Kreis, oder anders formuliert, Deine Aufruffrequenz wird durch 40.265 
geteilt.

Bei den Filtern muss man aufpassen, dass nichts aufschwingt, gerade bei 
Integerrechnung wichtig. Es gibt auch spezielle Filterstrukturen, die 
besser für die Erzeugung geeignet sind als die verwendete. Bei Nachfrage 
gern mehr.

Cheers
Detlef

von Elektroingenieur (Gast)


Lesenswert?

Das Aufschwingen lässt sich verhindern, indem man Verluste einbaut. Ein 
fach einige Promille der Werte versickern lassen. Das stabilsiert.

von Klausi (Gast)


Lesenswert?

Hallo,
schonmal vielen Dank an die Antworter.

Ich komme so langsam mit.
Ich habe den Generator angeschmissen in dem ich wie unten zu sehen das 
Array- Feld[1] mit einem Wert vorlade, und dann einfach sich selbst 
überlassen. Das funktioniert auch ganz fein. Ich wurde dann aber 
neugieriger und wollte die Frequenz einstellbar machen und möglichst 
einen "perfekten" Sinus erhalten.

Der erste Versuch war, eine Verzögerung einzubauen, was aber nicht sehr 
schön ist. Hier der Code dazu, bitte seht über die unsinnige 
Parameter-Übergabe oder Nicht-Übergabe hinweg. Für Tipps und Vorschläge 
wie ich oben beschriebenes gut implementieren könnte, wäre ich sehr 
froh:

int IIR_sin_gen(int delay){
 int i;
 //for (i=0; i<delay; i++) asm("nop;");// Verzögerung

 int output;
 static int A = 0x7e66;

 //pi/40->0x7e66
 // pi/1meg->0xffffac1b ;
 // A=(1.975/2 * 32768)

 static int y[3] = { 0, 0x120902de, 0};

 //y1->pi/40->0x120902de
 //y1=(0.1409*32768)

 y[0] = (((A*y[1])>>15) + ((A*y[1])>>15)) - y[2];
 y[2] = y[1];   /* y2 <-- y1 */
 y[1] = y[0];    /* y1 <-- y0 */
 output = y[0];
 return ( output);

} //IIR_sin_gen()

von Leo (Gast)


Lesenswert?

Hallo Klausi,

ich finde deinen Beitrag sehr interessant.Wie kommst du auf diese 
Gleichung? --> y[n]= A*y[n-1] + B*y[n-2] + x[n]

Und warum 2. Ordnung ?

Man könnte doch eigentlich einen IIR Filter 1. Ordnung realisieren, der 
nur einen Pol auf dem Einheitskreis hat ? Dies führt auch dazu, dass das 
System dann schwingt, oder ?

von detlef_a (Gast)


Lesenswert?

1-----
>>Der erste Versuch war, eine Verzögerung einzubauen, was aber nicht sehr
>>schön ist.

Du baust keine weitere Verzögerung ein, du machst damit kein Filter 3. 
Ordnung, das y[0] benutzt Du nur als Zwischenspeicher ohne es in der 
Rückkopplung zu verwenden, Dein Filter ist nach wie vor 2.Ordnung und 
das reicht auch, 1.Ordnung reicht zu Schwingen nicht. Vergiss das y[2], 
y[0] und y[1] reicht.

2----
y initialisierst Du mit 0x120902de, dann macht Dir ((A*y[1])>>15)) nen 
Overflow bei 32 Bit int Rechnung. Bei 32 Bit Rechnung darf A max. 16Bit 
sein (ist es) und die y auch (die sind größer).

3----
>>>Ich habe den Generator angeschmissen in dem ich wie unten zu sehen das
>>>Array- Feld[1] mit einem Wert vorlade, und dann einfach sich selbst
>>>>überlassen.

Das ist richtig so: y 'vorladen', das bestimmt Dir Amplitude und Phase 
des Generators und hat mit der Frequenz nix zum tun. Die Frequenz wird 
einzig und allein durch A bestimmt. Dein Filter ist stabil und richtig, 
wenn Du (in Deinem Fall) nach 40-Umläufen wieder die Initialwert in y[0] 
und y[1] stehen hast. Falls nein muß Du das erzwingen oder genauer 
rechnen (64 Bit) oder mit der Genauigkeit des A fummeln.

4----
Das geht besser: Die Größen y[0] und y[1] betrachtest Du als die 
Komponenten eines Vektors, den Du um 2*pi/40 drehst.
So:
y0alt=y[0];
y1alt=y[1];
y[0]=cos(2*pi/40)*y0alt-sin(2*pi/40)*y1alt;
y[1]=sin(2*pi/40)*y0alt+cos(2*pi/40)*y1alt;
Alles natürlich auf Integerrechnung umlügen.
Das kostet mehr Multiplikationen (ich glaube 2 statt einer wenn man 
Ergebnisse clever zwischenspeichert) ist aber numerisch stabiler. Kann 
man auch besser 'vorladen', nämlich direkt mit der gewünschten 
Amplitude.
http://de.wikipedia.org/wiki/Rotationsmatrix


math rulez!
Cheers
Detlef

von Detlef _. (detlef_a)


Lesenswert?

Ich brauchte Spaß und habe die Rotation in Integerrechnung gebastelt, 
Sourcecode angehängt.
Das implementiert einen Quadraturgenerator (Sinus sind die y[0], Cosinus 
sind die y[1] ), der nach 40 Iterationen einmal GENAU auf seinen 
Ausgangspunkt zurückkommt). y[0] habe ich immer mit 0 initialisiert, 
Initialisierung von y[1] habe ich möglichst groß machen wollen. Die 
defines der SI und CO habe ich mit Korrekturwerten versehen, +1 bzw. +2, 
sodass der Initialisierungswert des y[1] möglichst groß wurde. Dann war 
22370 das Maximum, das ich gefunden habe. Das ist auch die einzige 
Lösung, die ich gefunden habe, die mit Initialisierung mit einer 2er 
Potenz funktioniert, nämlich 1024.
Das ist eine anschauliche Implementierung, das gleiche läßt sich 
effektiver coden.

Cheers
Detlef

int16_t y[2];
int16_t n,k,m;
int16_t yalt[2];

#define SI ((int16_t)(32767*sin(2*3.141592654/40))+1)
#define CO ((int16_t)(32767*cos(2*3.141592654/40))+2)

y[0]=0;
y[1]=22370;
//y[1]=1024;
for(k=0;k<40;k++){
 yalt[0]=y[0];
 yalt[1]=y[1];
 y[0]= (int16_t)(((int32_t)yalt[0]*CO-(int32_t)yalt[1]*SI)>>15);
 y[1]= (int16_t)(((int32_t)yalt[0]*SI+(int32_t)yalt[1]*CO)>>15);
 printf("k %d y[0] %d y[1] %d\n",k,y[0],y[1]);
}
return;

von Klausi (Gast)


Lesenswert?

Wow, vielen Dank an Detlef_a, ich werde das erst morgen ausprobieren und 
nachvollziehen können. Bis Dahin viele Grüße

von Dieter (Gast)


Lesenswert?

Hallo zusammen,

sehr interessanter Beitrag ich habe im Moment die Aufgabe ein 
Synthesizer auf zubauen vielleicht auch mit der oben beschriebenen 
Methode. Nun habe ich jedoch eine Frage gibt es eine Möglichkeit um die 
Frequenz zu ändern damit ich ein Sweep habe?

von Sascha (Gast)


Lesenswert?

Was sehr einfaches ist hier:
http://www.musicdsp.org/archive.php?classid=1#10

Ist vom Prinzip her ein State-Variable-Filter.

Der Coeffizient 'a' muss ja nicht pro Sample aktualisiert werden, wenn 
man blockweise arbeitet. Z.B. alle 16 oder 32 Samples, das sollte für 
einen Sweep reichen.

von BiBi (Gast)


Lesenswert?

Wenn es möglichst kontinuierlich sein sollte, wird man schon pro Sample 
arbeiten müssen, oder?

von Sascha (Gast)


Lesenswert?

Bei 44k SR wäre das 'ne halbe ms pro Update. Das hört kein Mensch. 
Selbst für 'zap'-Sounds (von oben nach unten) ist das noch schnell 
genug. Die meisten Software-Synths machen das nur alle 16 - 64 Samples.

Eine andere Lösung wären halt normale Biquads. Aber die sollten für 
schnelle Sweeps in Direct Form I sein, nicht II, andernfalls werden die 
bei schnellen Koeffizenten-Updates schnell instabil.

Klar, das Update kann auch pro Sample sein. 1x Sinus ist ja auf 
aktuellen Systemen nicht mehr so teuer (LUT würde ich wg. 
Cache/Speicherzugriff nur noch selten machen).
Bei den Chamberlins sollte man aber beachten, dass man mindestens 4x 
oversamplen sollte, um bei 44k den gesamten hörbaren Bereich 
durchsweepen zu können.

von J. S. (grooc)


Lesenswert?

Hallo, Kann mir einer nochmal genauer erklären wie man auf die 
Startwerte kommt. Mir ist aufgefallen das es nicht mit allen Werten 
funktioniert, aber ich bin mir sicher, dass es eine Abhängigkeit gibt, 
die das Problem Mathematisch löst oder?

Habe es jetzt gefunden!
hier steht es:  http://www.musicdsp.org/archive.php?classid=1#10

Init:
w = omega = f/fa*pi
ip = initial phase e.g. pi/4
y1=sin(ip-w)
y2=sin(ip-2*w)

: Bearbeitet durch User
von J. S. (engineer) Benutzerseite


Lesenswert?

Man kann sich das am Einfachsten so vorstellen, dass es ein idealer 
Oszillator ist, der irgendeine Form von Anfangsenergie benötigt, die 
dann das Filter selbständig schwingen lässt. Die Methodik krankt aber 
nach meinen Erfahrungen an der Akkumulation der Fehler. Wenn man kein 
natürliches Verhalten eines Oszillators braucht, ist DDS letztlich 
einfacher und genauer.

von J. S. (grooc)


Lesenswert?

Noch eine Frage :) jetzt habe ich die tolle Formel gefunden, aber wie 
stelle ich damit die Amplitude ein? oder wähle ich die Phase einfach so 
dass ich als Startwerte möglichst große Werte habe?

von Sascha (Gast)


Lesenswert?

J. S. schrieb:
> aber wie
> stelle ich damit die Amplitude ein? oder wähle ich die Phase einfach so
> dass ich als Startwerte möglichst große Werte habe?

Steht doch dort in den comments:
"FYI you can set s[0] to whatever amplitude of sinewave you desire. With 
0.5, you will get +/- 0.5"

von T.U.Darmstadt (Gast)


Lesenswert?

Warum multiplisierst Du nicht einfach den Ausgangswert mit dem 
benötigten Faktor? Den Sinusgenerator würde ich lieber auf maximaler 
Auflösung rechnen- und schwingen lassen

von chris_ (Gast)


Lesenswert?

Hier hätte ich nocn einen einfachen Sinusgenerator:
http://www.hobby-roboter.de/forum/viewtopic.php?f=5&t=149

von T.U.Darmstadt (Gast)


Lesenswert?

Der entspricht aber weitgehend dem im Beitrag weiter oben, nur das 
einige Parameter weggefallen sind, weshalb man umsoeher skalieren muss.

von chris_ (Gast)


Lesenswert?

Du meinst den hier:
Beitrag "Re: iir Filter als Sinusgenerator"
Da sind zwei Multiplikationen drinn. In meiner Version keine. Auf einem 
AVR ist da ein exorbitanter Geschwindigkeitsunterschied.

von J. S. (engineer) Benutzerseite


Lesenswert?

Beim Vergleich der ganzen Lösungen muss man teilweise etwas genauer 
hinsehen: Nicht alles, was sich direkt mathematisch formulieren lässt, 
ist hinterher auch schnell, zumal ein händischer Sinus noch etwas 
anderes ist, als ein vollständiger Sinusoszillator. In Sachen 
Genauigkeit gibt es da auch Abweichungen und natürlich macht es noch 
einen Unterschied, ob man das auf einem mini-Controller, einem FP-DSP 
oder einen FPGA realisiert.

Hier wäre eine manuelle Implementierung der verketteten Integratoren von 
oben fürs FPGA der einen realistischen Einschwingvorgang abarbeiten 
kann, mit Eingriffen für Frequenz, Beschleunigung, Dämpfung und 
Dithering für hochgenaues Ausschwingen auf Null.

http://96khz.org/doc/vasynthesisvhdl

von Tim  . (cpldcpu)


Lesenswert?

chris_ schrieb:
> Du meinst den hier:
> Beitrag "Re: iir Filter als Sinusgenerator"
> Da sind zwei Multiplikationen drinn. In meiner Version keine. Auf einem
> AVR ist da ein exorbitanter Geschwindigkeitsunterschied.

Auch ein Shift ist eine Multiplikation...

von Matthias (Gast)


Lesenswert?

...aber bei µC mit einem Barrelshifter nur 1 Cycle zum Abarbeiten (zB 
dsPICs haben sowas eingebaut).

von chris_ (Gast)


Lesenswert?

... und bei einem ARM-Prozessor der Barrelshifter direkt hinger der ALU 
hängt und keine Rechenzeit frist.

http://www-mdp.eng.cam.ac.uk/web/library/enginfo/mdp_micro/lecture4/lecture4-3-2.html

von experte (Gast)


Lesenswert?

Im FPGA würde der auch nichts verbrauchen.

Bitte melde dich an um einen Beitrag zu schreiben. Anmeldung ist kostenlos und dauert nur eine Minute.
Bestehender Account
Schon ein Account bei Google/GoogleMail? Keine Anmeldung erforderlich!
Mit Google-Account einloggen
Noch kein Account? Hier anmelden.