Post by Johannes BauerEs geht nicht um irgendwie physikalisch korrektes rosa/braunes Rauschen,
sondern darum, dass das -- klingt blöd, weiß ich -- "sanft" klingt und
eben in Software "abschwächbar" ist. Konkret verändere ich also die
Verteilung Energie/Frequenz in Software.
Das ist unheimlich simpel: Ich habe einen 8-Bit Akkumulator, der jeweils
auf PORTA ausgegeben wird. In jeder Iteration wird der durch
accu += (prng() / divisor)
Der Divisor hier ist eigentlich Grütze. Er reduziert nur die Zahl der
Zufallsbits, die von prng kommen. Sehr deutlich wird das bei
Zweierpotenzen. Diese entsprechen einem Shift-Operator, der einfach die
unteren Bits entsorgt.
Post by Johannes BauerVoller Code: https://github.com/johndoe31415/tiny26noise/blob/master/main.c
Verändert. Je größer "divisor", desto "sanfter" das Rauschen.
Prinzipiell funktioniert das auch prima, mit einer Einschränkung: In der
Praxis erhalte ich mit größeren Divisor-Werten auch immer ein lästiges
knacksen.
Überläufe?
Post by Johannes BauerDa das immer an denselben Stellen auftritt in demselben Pattern (ist ja
ein 100% deterministischer PRNG) nehme ich an, dass das an der Division
liegt, die auf dem AVR in Software emuliert wird und daher nicht
Laufzeitkonstant ist.
Das würde mich wundern, wenn das der wesentliche Punkt ist, zumal der
Divisor ja auch noch konstant ist.
Berechne doch mal die Folge auf dem PC vor und schau dir das Ergebnis
mit Audacity an.
Post by Johannes BauerJetzt habe ich überlegt, wie ich das besser machen kann. Klar, kann nur
Zweierpotenzen als Teiler zulassen, aber dann ist der Unterschied
zwischen zwei Levels recht stark. Fällt jemandem eine gute Idee ein, wie
man das ohne Hardwaredivision, nur mit Integerarithmetik und in
konstanter Laufzeit hinkriegt? Mir fehlt gerade eine gute Idee dafür.
Ich würde komplett anders vorgehen.
prng() liefert die Pseudo-Zufallsbits, also Weißes Rauschen. Dieses
würde ich durch ein IIR-Filter jagen und damit passend maskieren. Im
einfachsten Fall ist das ein Integrator (Addition), dann gibt es erst
mal Rosa Rauschen. Du musst nur aufpassen, dass Du keine Überläufe
bekommst, sonst ist das Spektrum komplett im Eimer.
Das ist deiner Methode bis dahin gar nicht so unähnlich (die Addition),
bis auf die Sache mit den Überläufen eben.
Jetzt kann man natürlich den Erwartungswert von prng() bei jeder
Addition abziehen und somit eine Eskalation des Akkumulators
abschwächen, verhindern wird es sie aber nicht, denn Rosa Rauschen hat
nun einmal bei der Frequenz 0 eine Singularität (unendliche Energie).
Technisch wäre es zudem eine Null-Operation, da es nur das oberste Bit
von prng toggeln würde, das ohnehin zufällig ist. Also bringt nichts.
Das bringt mich zum zweiten Punkt: in jedem Fall braucht der Akkumulator
mehr Bits (also faktisch 16), von denen nur die obersten an den ADC
gehen. Jetzt hat man Platz für weitere Kalkulationen.
Jetzt kommt der eigentliche Clou: ein geeignetes Filter würde die nach
obigem Muster kompensierten Akkumulatorwerte zwar ebenfalls
tiefpassfiltern, aber es hat eben auch eine Resonanzfrequenz unterhalb
der nicht mehr verstärkt oder besser noch Hochpass filtert. Damit ist
man den lästigen DC-Anteil los und ebenso sind unnötige Rumpelgeräusche
weg, die letztlich DAC-Headroom kosten, von denen der 8-Bitter ohnehin
viel zu wenig hat.
Man braucht also am besten einen Bandpass mit einer sehr schlechten
Güte, vorzugsweise aperiodisch, also Q=0,7. Die Resonanzfrequenz des
Filters ist grob die gewünschte untere Grenzfrequenz des Rosa Rauschens.
Technisch ist ein IIR Filter ein bisschen +,-,*. Division braucht man
dafür gar nicht. Man nennt das auch Biquad, weil es sich mit einer
biquadratischen Gleichung modellieren lässt.
y[n] = a0 x[n] + a1 x[n-1] + a2 x[n-2] + b1 y[n-1] + b2 y[n-2]
y sind die Ausgabewerte, x die Eingabewerte aus prng. [n-1] meint den
letzten Ein- bzw. Ausgabewert, [n-2] den vorletzten.
Die Koeffizienten a0 .. b2 bekommt man z.B. hier her:
https://www.earlevel.com/main/2013/10/13/biquad-calculator-v2/
Dort "Bandpass" wählen und Güte 0,7 (default).
Die Koeffizienten hängen nur von dem Verhältnis untere Grenzfrequenz zu
Samplingrate bzw. Schleifenfrequenz ab.
Bei allen hier relevanten Filter ist a1 = 0 - der Term entfällt also und
man braucht nur noch den vorletzten Eingabewert - Und a2 = -a0, zudem
ist b1 immer negativ. Damit reduziert sich die Formel zu
y[n] = a0 x[n] - a0 x[n-2] - (-b1) y[n-1] + b2 y[n-2]
Da x[n] ohnehin zufällig ist, ist es auch x[n] - x[n-2]. Man kann sich
also die Nummer mit x[n-2] komplett schenken.
y[n] = a0 x[n] - (-b1) y[n-1] + b2 y[n-2]
Die drei y-Akkumulatoren müssen wie schon erwähnt 16-bittig ausgeführt
werden. Für x genügen 8 Bit, mehr kommt von prng ohnehin nicht.
Die Koeffizienten liegen aufgrund des großen Verhältnisses von Unterer
zu Samplingrate alle nahe bei ganzen Zahlen. Daher bietet es sich an,
die Koeffizienten ebenfalls nur 8-Bittig auszuführen. Gespeichert wird
dann nur die Abweichung zur nächsten ganzen Zahl. Damit reduzieren sich
die Multiplikationen für die auf 8 Bit mal 16 Bit.
Beispiel: b1 = -1,98; per Definition sind die oberen 8 Bit des
y-Akkumulator die Vorkommastellen, die später an den DAC gehen, und die
unteren 8 Bit die Nachkommastellen.
b1 * y[n-1] = b1' * y[n-1] - (y[n-1] << 1) mit
b1' = b1 + 2 = 0,02 entspricht 0x05
Das ergibt natürlich eine sehr grobe Quantisierung (also die 0x05).
Anders formuliert, auch 16 Bit reichen nicht. Das Filter wäre instabil.
Man braucht mehr.
Man kann sich aber leicht mehr Nachkommastellen herauskitzeln, indem man
mit Shift-Operatoren arbeitet:
b1 * y[n-1] = ((b1' * y[n-1]) >> 5) - (y[n-1] << 1) mit
b1' = 32 * 2 - b1 = 0,64 entspricht 0xA4
Analog geht man für die anderen Terme vor.
Bleibt noch die Sache mit der Pegelregelung. Dazu nimmt man sich den
a0-Koeffizienten vor. Also:
a0 * x[n] = (a0' * x[n]) >> S
Man beachte die Shift-Richtung. Per Definition sind x[n]
Vorkommastellen, also die oberen 8 Bit. Diese werden multipliziert mit
den Nachkommastellen von a0. Heißt, die oberen 8 Bit des 16-bittigen
Multiplikationsergebnisses sind am Ende Vorkommastellen, die unteren
Nachkomma, exakt passend zu dem y-Akkumulator.
Allerdings wird man schnell feststellen, dass der y-Akkumulator auf
diese Weise überläuft. Daher muss das Signal noch zusätzlich
abgeschwächt werden. Das macht man mit S. Zudem muss man sich auch hier
für a0 zusätzliche Bits erarbeiten, analog zu b1 oben. Daher wird S noch
größer. Vielleicht wählt man S = 8, was nichts anderes bedeutet, als die
oberen 8 Bit des Multiplikationsergebnisses werden nur zu den
Nachkommastellen des y-Akkumulators addiert. Damit muss man nur noch bei
Übertrag an die oberen 8 Bit von y ran.
Die eigentliche Pegelanpassung mach man jetzt über S (in 6dB-Schritten)
und dem a0'-Koeffizienten für's feine. Dieser hat keinen Einfluss auf
den Frequenzgang.
Am Ende muss man ausprobieren, wie weit hoch man sich traut. Es ist
dabei essentiell, dass der y-Akkumulator niemals überläuft.
Alternativ könnte es auch klappen, wenn S = 8 und a0' = 1.0 wählt, also
gar keine Multiplikation mit a0. Damit bekommt man die weiter reduzierte
Formel:
y[n] = (x[n] >> 8)
+ ((b1' * y[n-1]) >> 5) - (y[n-1] << 1)
- ((b2' * y[n-2]) >> 5) + y[n-2]
wobei >> 8 jetzt nur der Addition zu den unteren 8 Bit von y entspricht.
Man braucht jetzt genau 2 8*16 Bit Multiplikationen pro Ausgabewert und
keine Division.
Alle bisher genannten Kalkulationen beziehen sich auf Signed-Integer
Werte. Technisch hat das bis hier hin keine Auswirkung, aber bei der
Übergabe an den DAC muss man das beachten. hier genügt es nicht, die
oberen 8 Bit des y-Akkumulators zu nehmen, man muss zudem das oberste
Bit toggeln (entspricht +128).
Kleiner Hinweis noch: Die Anzahl der signifikanten Bits bei der
Kalkulation der b-Koeffizienten sind eklatant wichtig für Genauigkeit
und Stabilität des Filters. Die beiden b-Terme liefern stets ähnlich
große Zahlen, die sich gegenseitig weitgehend Auslöschen. Stimmt das
Verhältnis der beiden durch Quantisierungsfehler nicht, so wächst y
schnell mal exponentiell über alle Grenzen.
Während man peinlich darauf achten muss dass y[] iterativ nicht
überläuft, sind temporäre Überläufe während der Kalkulation eines y[n],
beispielsweise durch y[n-1] << 1 hingegen unkritisch. Die kompensieren
sich automatisch.
Marcel