Forum: Mikrocontroller und Digitale Elektronik Samplerate conversion - rationaler Faktor


von Vincent H. (vinci)


Lesenswert?

Grüß euch

Ich bin auf der Suche nach effizienten Algorithmen zur Konvertierung von 
Audio-Sampleraten um rationale Faktoren. Diese sollen auf einem STM32L4 
für 16bit/44kHz Samples in Echtzeit laufen, vorzugsweise auf mehr als 
einem Dutzend Kanälen...

Meine Recherche hat ergeben, dass hierfür wohl meist kaiser windows in 
Form von FIR Filtern genutzt werden. Siehe z.B. hier:
https://ccrma.stanford.edu/~jos/resample/resample.pdf
Bzw. als Library hier:
https://github.com/depp/libfresample

Leider bin ich nicht sicher wie realistisch diese Methode für meine 
Anwendungszwecke ist. Grad unter dem Link der Library findet sich ein 
Beispiel in dem erwähnt wird, dass ein 1.6 GHz Intel Atom für eine 
16bit/48->44.1 Wandlung in "mittlerer Qualität nur" etwa 197x schneller 
als Echtzeit rechnet. Jetzt kann ich als Laie aus dem Stehgreif schwer 
sagen was "mittlere Qualität" in Bezug auf die Taps-Länge heißt, aber 
wenn ein 1.6GHz Atom "nur" 197x schneller als Echtzeit ist, mach ich mir 
bei meinem 80MHz L4 schon irgendwie Sorgen...

Leider ist die von ST Angebote Library zur Umwandlung von 44.1->48kHz 
nicht Open Source... Bieten andere Hersteller da irgendwas an wo man 
reinschaun kann?

mfg
Vincent

/edit
Pardon, hät wohl besser in "DSP" gepasst

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


Lesenswert?

Der saubere Weg in dem Fall ist eine 441 zu 480 Ratenumsetzung, bzw 
umgekehrt mit einem CIC-Filter. Für höchste Qualität vor und nach dem 
CIC ein FIR-Filter mit Eckfrequenz von z.B. -6dB@20kHz. Wenn man sparen 
will, nimmt man die beiden Frequenzen asymmetrisch, d.h die eine enger 
am Band, die andere weiter und zudem ganzzahlig zu der Abtastfrequenz, 
also im anderen Fall 44,1 / 2 auf 48kHz und 48kHz/3 bei 44kHz. Das ist 
aber für PCs zu aufwändig.

Mit einigen Bauchschmerzen geht daher auch 11/12 und Filterung über 3 
samples, also zwei FIR-Filter mit 3x11 oder 3x12 TAPs und der 
entsprechenden Frequenz. Tiefpass-Filter-Koeffizienten aus Matlab.

von Jobst M. (jobstens-de)


Lesenswert?

Die letzten zwei Samples des 44100-Signals werden Gespeichert.

In einer Tabelle liegen die Integralwerte einer vollen, in 160 Stücke 
unterteilten, invertierten Cosinuswelle auf 1 normiert. (Was eine halbe, 
invertierte Cosinuswelle von 0 bis 1 ist)
(Oder die Integralwerte der Filterkoeffizienten, die man benutzen 
möchte. Die virtuelle Samplingfrequenz liegt bei 7,056MHz.)

Von einem Zeiger auf die Tabelle wird nach jedem 48000-Sample 13 
abgezogen und ein neues 44100-Sample geladen. Ist der Zeiger negativ, 
werden 160 Addiert und KEIN neues 44100-Sample geladen.

Der neuere der beiden 44100-Samples wird mit dem Wert aus der Tabelle 
multipliziert, der ältere der beiden 44100-Samples wird mit (1 - Wert 
aus der Tabelle) multipliziert. Die beiden Ergebnisse werden Addiert und 
sind das fertige 48000-Sample.

Das Ganze lässt sich in 2x 16 Bit in unter 100 Takten in z.B. einem 
AT-Mega durchrechnen. Mit 7,056MHz takten, einen Timer auf 147 Takte, 
den zweiten auf 160 Takte, dann hat man auch den Input/Output-Takt.
Auf Prozessoren, welche zwei 16-Bit Zahlen am Stück multiplizieren 
können geht es flotter. Mit höherem Takt auch, so dass ich auf Deiner 
80MHz ARM CPU für ein Duzend Kanäle gar keine Probleme sehe.


Gruß

Jobst

von Vincent H. (vinci)


Lesenswert?

Jürgen S. schrieb:
> Der saubere Weg in dem Fall ist eine 441 zu 480 Ratenumsetzung, bzw
> umgekehrt mit einem CIC-Filter. Für höchste Qualität vor und nach dem
> CIC ein FIR-Filter mit Eckfrequenz von z.B. -6dB@20kHz. Wenn man sparen
> will, nimmt man die beiden Frequenzen asymmetrisch, d.h die eine enger
> am Band, die andere weiter und zudem ganzzahlig zu der Abtastfrequenz,
> also im anderen Fall 44,1 / 2 auf 48kHz und 48kHz/3 bei 44kHz. Das ist
> aber für PCs zu aufwändig.
>
> Mit einigen Bauchschmerzen geht daher auch 11/12 und Filterung über 3
> samples, also zwei FIR-Filter mit 3x11 oder 3x12 TAPs und der
> entsprechenden Frequenz. Tiefpass-Filter-Koeffizienten aus Matlab.


Hast du vielleicht ein paar Literaturtips zur CIC Umsetzung für mich? Da 
gibts vor allem in Bezug auf embedded noch nicht sonderlich viel find 
ich. Brauch ich bei der Realisierung eines CIC einen Buffer in der 
vollen Größe? Sprich eigentlichte Buffergröße zum Abspielen * 
Upsample-Faktor (z.B. 441)?



Jobst M. schrieb:
> Die letzten zwei Samples des 44100-Signals werden Gespeichert.
>
> In einer Tabelle liegen die Integralwerte einer vollen, in 160 Stücke
> unterteilten, invertierten Cosinuswelle auf 1 normiert. (Was eine halbe,
> invertierte Cosinuswelle von 0 bis 1 ist)
> (Oder die Integralwerte der Filterkoeffizienten, die man benutzen
> möchte. Die virtuelle Samplingfrequenz liegt bei 7,056MHz.)
>
> Von einem Zeiger auf die Tabelle wird nach jedem 48000-Sample 13
> abgezogen und ein neues 44100-Sample geladen. Ist der Zeiger negativ,
> werden 160 Addiert und KEIN neues 44100-Sample geladen.
>
> Der neuere der beiden 44100-Samples wird mit dem Wert aus der Tabelle
> multipliziert, der ältere der beiden 44100-Samples wird mit (1 - Wert
> aus der Tabelle) multipliziert. Die beiden Ergebnisse werden Addiert und
> sind das fertige 48000-Sample.
>
> Das Ganze lässt sich in 2x 16 Bit in unter 100 Takten in z.B. einem
> AT-Mega durchrechnen. Mit 7,056MHz takten, einen Timer auf 147 Takte,
> den zweiten auf 160 Takte, dann hat man auch den Input/Output-Takt.
> Auf Prozessoren, welche zwei 16-Bit Zahlen am Stück multiplizieren
> können geht es flotter. Mit höherem Takt auch, so dass ich auf Deiner
> 80MHz ARM CPU für ein Duzend Kanäle gar keine Probleme sehe.
>
>
> Gruß
>
> Jobst

Interessante Methode, aber wohl nur auf diesen Spezialfall anwendbar? 
Ich benötige leider eine Methode, die fließend auf x-beliebige Fälle 
anwendbar ist und sich nicht auf 44.1->48 beschränkt. Grad dies 
erschwert die Sache zusätzlich, da man wohl für jeden Fall neue 
Filterkoeffizienten bräuchte, ganz egal ob nun FIR/IIR oder sonst was?

von Jobst M. (jobstens-de)


Lesenswert?

Vincent H. schrieb:
> Interessante Methode, aber wohl nur auf diesen Spezialfall anwendbar?
> Ich benötige leider eine Methode, die fließend auf x-beliebige Fälle
> anwendbar ist und sich nicht auf 44.1->48 beschränkt. Grad dies
> erschwert die Sache zusätzlich, da man wohl für jeden Fall neue
> Filterkoeffizienten bräuchte, ganz egal ob nun FIR/IIR oder sonst was?

Fließend wird damit schwierig. Dynamisch ist möglich, indem man die 
Restfaktoren, welche nicht gemeinsam sind, errechnet und damit die 
Tabelle neu berechnet.


Für wirklich fließende (also FS ändert sich stetig) wird Dir nur 
bleiben, ein langes Register mit Werten mit hohem Takt anzulegen und bei 
Bedarf (Ausgangssample wird benötigt) eine Berechnung durchzuführen. 
Quasi ein asyncrones FIR-Filter.


Gruß

Jobst

von Dergute W. (derguteweka)


Lesenswert?

Moin,

Zu den CIC-Filtern gibts hier was:
http://www.dspguru.com/dsp/tutorials/cic-filter-introduction

Sind aber imho fuer eine CPU nicht gut geeignet, weil bei "heutigen" 
CPUs Multiplikationen nicht viel schmerzhafter als Additionen sind, 
dafuer aber eine vergroesserung der Bitbreite mehr Schmerzen macht.
CIC sind "schoen" in FPGAs.

Der STM32 kann doch "DSPisch" - also wird's auf eine 
Polyphasenfilterbank bzw. ein einzelnes FIR, das sich seine 
Koeffizienten jeweils aus einer Tabelle holt, rauslaufen. Die Anzahl 
Taps des FIR bestimmt dann benoetigte Rechenleistung und Qualitaet der 
Konversion.

Das Kochrezept von Jobst M laeuft auf sowas raus - mit einem 2-Tap-FIR 
entsprechend Rechenleistungsschonend.

Gruss
WK

von Bla (Gast)


Lesenswert?

Vincent H. schrieb:
> Interessante Methode, aber wohl nur auf diesen Spezialfall anwendbar?
> Ich benötige leider eine Methode, die fließend auf x-beliebige Fälle
> anwendbar ist und sich nicht auf 44.1->48 beschränkt. Grad dies
> erschwert die Sache zusätzlich, da man wohl für jeden Fall neue
> Filterkoeffizienten bräuchte, ganz egal ob nun FIR/IIR oder sonst was?

Das geht im Grunde auch so. Basierend auf seinem Beitrag hab ich gestern 
damit verbracht, mir eine komplett unabhängige ASRC zu basteln. Die 
Performance ist (für den geringen Rechenaufwand) echt gut! Ok, kein 
sinc, aber dafür auch deutlich geringere Latenz!

Du brauchst für eine ASRC nach diesem Prinzip:
-> einen möglichst genauen Kosinus, einmal pro Sample. Also eher keine 
LUT sondern irgendwie berechnen.
-> zwei Speicherworte für das letzte und vorletzte x und x_1
-> einen "hinreichend genauen" Zähler c (ich würde 16 bit nehmen)
-> das Verhältnis df = f_in/f_out (Capture-Modul oder Frequenzzähler)

Für jede Ausgangssample:
gewicht = (cos(c)+1)/2;
out = x * gewicht + x_1 * (1-gewicht);
c = c + df;

Der Trick liegt darin, schlau zu skalieren.
c und cos sollten so skaliert sein, dass cos(c=0) = 1 ist und 
cos(c=65535) = -1
df wird mit 65536 multipliziert.

Das einzig kompliziertere ist der Kosinus, ansonsten läuft das relativ 
präzise ohne Modulo, Vergleiche, Sprünge.

Für den Kosinus gibts auch einen Trick: Wenn man den Kosinus so 
berechnet wie hier: 
http://www.earlevel.com/main/2003/03/02/the-digital-state-variable-filter/, 
d.h. als Phasenschieber-Oszillator, kann man statt cos(c) auch mit dem 
dortigen code, und f = df, berechnen. Da gibts nur das Problem, dass wir 
den Kosinus nur von 0 bis pi wollen, wenn wir das aber mit dem obigen 
Oszillator berechnen, läuft unser Kosinus von 0 bis 2*pi statt 0-pi. Wir 
müssen also nach der Hälfte eine Kosinus-Zyklus um eine halbe Phase nach 
hinten springen - oder, dank Periodizität des Kosinus, ihn einfach 
invertieren, d.h statt x * gewicht + x_1 * (1-gewicht) eben x * 
(1-gewicht) + x_1 * gewicht rechnen.

Wir müssen also nach jedem halben Zyklus den Kosinus invertieren. Dazu 
schauen wir uns einfach das Vorzeichen des Sinus an, das wechselt schön 
alle halbe Zyklus phasenkorrekt.

Folgender Code sollte daher theoretisch simpel und schön eine ASRC 
berechnen:

(Ungetesteter Code, bei der signed-unsigned-Multiplikation muss evtl. 
sauberer gecastet werden)
1
int16 sinz = 0;
2
int16 cosz = 32767;
3
int16 in, in_prev;
4
int16 df = 0;
5
uint16 in_clocks = 0;
6
uint16 out_clocks = 0;
7
8
update_input(int16 new_in){
9
  in_clocks++;
10
  if(in_clocks==0){
11
    df = (int16)out_clocks;
12
    out_clock = 0;
13
  }
14
  in_prev = in;
15
  in = new_in;
16
}
17
18
int16 update_output(){
19
  out_clocks++;
20
  sinz += df * cosz;
21
  cosz -= df * sinz;
22
  uint16 gewicht = (uint16)cosz + 0x8000;
23
  if(sinz>=0){
24
    return in * gewicht + in_prev * -gewicht;
25
  }else{
26
    return in * -gewicht + in_prev * gewicht;
27
  }

Wenns zu wirr ist, ich will mir das die Tage nochmal anschauen, weil ich 
das Konzept spannend finde!

df wird noch relativ langsam geupdated (alle 65ksamples, d.h. alle 1,5 
Sekunden). Da kann man über ne Divisionsroutine nachdenken, die df = 
out_clocks/in_clocks * 65636 ausrechnet, oder df alle 8192 samples 
updaten und einen Mittelwert der letzten df-Werte machen.

Unter Umständen ist df = in_clocks/out_clocks, dann muss der 
update-Block eben in update_output().

Weil ich grade so Spaß an Bitschubsen hab, hier noch ein möglicher 
modifizierter df-Code:
1
  if(!(in_clocks&0x1FFF)){    // alle 8192 Samples updaten
2
    int16 temp = (int16)out_clocks * 4;
3
    df = (df + temp)/2;       // df nicht so schnell ändern, sondern mitteln
4
    out_clock = 0;
5
  }

Mit dem Design sollten Frequenzunterschiede von 1/65536 auflösbar sein, 
d.h du hast ne Frequenzauflösung von 0.73Hz bei 48kHz.

von Bla (Gast)


Lesenswert?

1
   sinz += ((int32)df * cosz)>>16;
2
   cosz -= ((int32)df * sinz)>>16;

muss das natürlich heißen.

von Bla (Gast)


Lesenswert?

Ah, da merke ich, dass der Sprung vom Kosinus besser in der 
update_input() stattfinden sollte.

Alternativ kann man sich das Gespringe mit dem Sinus sparen, wenn man 
die zwei input samples als "Ringpuffer" betrachtet, d.h. abwechselnd 
updatet, und den Kosinus belässt wie er ist.

von J. S. (engineer) Benutzerseite


Lesenswert?

Dergute W. schrieb:
> Sind aber imho fuer eine CPU nicht gut geeignet, weil bei "heutigen"
> CPUs Multiplikationen nicht viel schmerzhafter als Additionen sind,
> dafuer aber eine vergroesserung der Bitbreite mehr Schmerzen macht.
> CIC sind "schoen" in FPGAs.

Das ist richtig, stimmt. Ich habe jetzt intuitiv an Hardware gedacht. 
Was man alternativ noch machen kann: Das Verfahren der einfachen 
Interpolation über 4 Punkte ausdehnen. Dazu gibt es verschiedene 
einfache Interpolationsverfahren.

von Vincent H. (vinci)


Lesenswert?

Vielen Dank für die vielen Inputs. Ich bin auf Grund der Komplexität des 
Themas noch nicht dazu gekommen viele Varianten durchzutesten... (schon 
gar nicht direkt auf der Hardware)

Außerdem fehlte mir bei den Filter-Zugängen teilweise etwas Background, 
da grad so Polyphase Sachen in verschiedenen Strukturen über "simple" 
IIR/FIR Filter hinaus gehn. Da war erstmal ein wenig einlesen angesagt 
(DSP using MATLAB und Multirate Filtering for DSP war da recht 
brauchbar).

Aus den beiden genannten Gründen hab ich zuerst mal MATLAB angeworfen 
und dort einige einfache Sachen ausprobiert. Unter anderem
- "zero order hold", sprich schlichtweg das vorrangegangene 
Original-Sample
- "nearest neighbor", sprich das nächstgelegende Original-Sample
- "lineare Interpolation"
- "cubische Interpolation"

und dann bereits wesentlich ausgefuchstere Sachen
- MATLABs implementiere re-sample Funktion (ein Upsampling gefolgt von 
einem FIR-Decimator angegebener Ordnung)
- ein selbstgebasteltes "Fractional Delay" Filter in Farrow-Structure 
(quasi auch ein FIR), dass wohl auch sehr hardware-freundlich sein soll, 
da sich unter anderem für verschiedene Delays die gleichen 
Filter-Koeffizienten nutzen lassen

Soweit so gut. Plottet man jetzt die Frequenzantwort jener Filter, bzw. 
gleich das Zeitsignal, dann sieht man natürlich wofür man den Aufwand 
treibt...

Und jetzt kommt das große aber.
Ich hörs nicht! Ich weiß ja dass ich kein audiophiler Mensch bin, aber 
auf meinem internen Billiglautsprecher am PC hör ich keinen bis kaum 
einen Unterschied zw. etwa "nearest neighbor" und FIR 4.Ordnung. Is das 
normal? Ich nehm mal an, dass durch die Hardware vom PC das Signal 
intern wieder "zurechtgesampelt" wird, egal was ich da für Sachen in 
MATLAB treib... aber dass man da so überhaupt nichts wahrnimmt?

Jetzt bin ich natürlich unschlüssig ob sich der Aufwand für mich 
überhaupt lohnen würde was anderes als linear oder übhaupt nicht zu 
interpolieren. Was ich anfangs vielleicht hätte erwähnen sollen ist, 
dass die Soundausgabe letztendlich auch nur auf sehr kleinen 
Lautsprechern stattfinden wird, wo man sich sowieso keine allzu gute 
Klangqualität erwartet...

Ich denke ich werde vorerst mal bei einem sehr simplen Algorithmus 
bleiben und erst bei Bedarf was "klügeres" implementieren. Ein großer 
Vorteil von "gar nix rechnen" ist übrigens, dass man sich keine Gedanken 
machen muss ob man jetzt extra Samples nur für die Berechnung in den 
aktuellen Buffer schieben muss... :D

von J. S. (engineer) Benutzerseite


Lesenswert?

... ganz klar: Die Unterschiede liegen ja bei den nichtharmonischen 
Verzerrungen und die billigen Lautsprecher dröhnen da mehr dazu, als 
sonst was. Nimm Dir mal gute Lautsprecher und einen sauberen Testton.

Und: Dank MP3 und oft zu lauter Musik haben viele Menschen schon Löcher 
im Gehörfrequenzgang. Daran liegt es zuweilen auch.

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.