Festkommarealisierung LMS (Least Mean Square)

[Mario A. (itse_me_mario)]

Hallo zusammen,

ich versuche gerade mein in Gleitkommaarithmetik entworfenes LMS-Filter 
in ein Festkomma LMS-Filter umzubauen. Leider habe ich im Internet noch 
keine Antwort auf mein Problem gefunden. Die Problemstellung ist in 
Zeichnung.png abgebildet. Vllt habt ihr auch noch Anmerkungen. Würdet 
ihr zus. eine Over-/Underflowerkennung einbauen? Würdet ihr long integer 
als Akkumulator verwenden?

Erst einmal habe ich CodeBlocks verwendet. Später soll das ganze noch 
auf einen DSP und Mikrocontroller sowie ein FPGA.

Das Gleitkomma-LMS-Filter ist im Anhang (siehe 
Referenzcode_LMS_Float_X1D1.txt). Die Funktion habe ich anhand von 
ERLE-Kurven überprüft - das passt. Zusätzlich habe ich das ganze auch in 
MATLAB implementiert.

Für ein Beta (Lernrate) von 0.01 habe ich die Ergebnisse für den 
Filterausgang y und den Fehler e angehängt (LMS_Data....txt). Damit man 
aber besser abschätzen kann, wo sich die Fehler im C-Programm befinden, 
habe ich die Signale x, d, e und y nochmal als .txt-Datei im Q15-Format 
angehängt.

Meine Daten d (DataD1Q15.txt) und x (DataX1Q15.txt) sind im Q1.15 
Datenformat. Dafür wurden die double Werte für d und x mit 2^15 
multipliziert. Die Ergebnisse sollen wieder als Q0.15 Zahlen (1 
Vorzeichen Bit, 15 Nachkommastellen) dargestellt werden.

Hier mein aktueller Stand, nach etlichen Änderungen... :P

#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>

#define N  32             // Anzahl adaptiver FIR-Filterkoeffizienten
#define NUM_ITERS  32000  // Anzahl Durchläufe

int16_t s16_beta = 327;   // Lernrate 0.01*2^15 = 327.68

int32_t s32_w[N] = {0};       
int16_t s16_w[N] = {0};       

int16_t s16_x[N] = {0};      

int32_t i = 0;
int32_t n = 0;

int16_t s16_d = 0;
int16_t s16_y = 0;
int16_t s16_e = 0;

int16_t s16_prod1 = 0;
int16_t s16_prod2 = 0;

int32_t s32_acc = 0;

int main()
{
    //*************************************************************
    FILE *fp_x;
    FILE *fp_d;
    FILE *fp_y;
    FILE *fp_e;
    //*************************************************************
    fp_x = fopen("DataX1Q15.txt","r");
    fp_d = fopen("DataD1Q15.txt","r");
    fp_y = fopen("LMS_Data_Y_Q15_beta_0_01.txt", "w");
    fp_e = fopen("LMS_Data_E_Q15_beta_0_01.txt", "w");
    //*************************************************************

    if( (fp_x == NULL) || (fp_d == NULL) || (fp_y == NULL) || (fp_e == NULL))
    {
        printf("One of the Files could not be opened.\n");
    }
    else
    {
        for(n=0; n<NUM_ITERS; n++)
        {
            //******************Read Data from TXT-File*************
            fscanf(fp_x,"%d\n", &s16_x[0]);
            fscanf(fp_d,"%d\n", &s16_d);
            //******************************************************

            s32_acc = 0; 

            for (i = 0; i < N; i++)
            {
                //hier könnte Overflow auftreten
                s32_acc = s32_acc + (int32_t)s16_w[i] * (int32_t)s16_x[i];
            }

            s16_y = (int16_t)(s32_acc >> 15);
            s16_e = s16_d - s16_y;

            for (i = N-1; i >= 0; i--)
            {
                s16_prod1 = (int16_t)((s16_x[i]  * s16_e   )>>15);
                s16_prod2 = (int16_t)((s16_prod1 * s16_beta)>>15);

                //Achtung hier könnte Overflow auftreten
                s16_w[i] = s16_w[i] + s16_prod2; 

                if(i!=0)
                {
                    s16_x[i] = s16_x[i-1]; 
                }
            }

            //************print values*************
            printf("%d Values: x=%d\td=%d\t Write Values \ty=%d\te=%d \t s32_acc=%d\n\n", n+1, s16_x[0], s16_d, s16_y, s16_e, s32_acc);
            //*************************************

            //******Write Data to TXT-File*********
            fprintf(fp_y, "%d\n", s16_y);
            fprintf(fp_e, "%d\n", s16_e);
            //*************************************
        }
    }
    fclose(fp_x);
    fclose(fp_d);
    fclose(fp_y);
    fclose(fp_e);

    while(1){}
}

Das Problem ist, dass die Werte für y sehr stark von d abweichen. Der 
Fehler e = d-y ist dadurch sehr groß und wird nicht kleiner mit der 
Zeit. Der Sollverlauf des Fehlers ist in 
LMS_Data_X1_D1_E_Float_Order31_beta_0_01.txt zu sehen. Also werden wohl 
die Filterkoeffizienten w[n] falsch bestimmt. Leider sehe ich nicht 
woran es liegt. Vermutlich sind die Werte für das rechtsshiften >>15 zu 
klein und werden dadurch 0 bzw. zu klein. Weiß jemand wie man das am 
besten angeht?

Später soll auch noch das NLMS Filter als Festkommavariante 
implementiert werden.

Über eure Hilfe oder/und Tipps für gute Literatur würde ich mich sehr 
freuen! Vielen Dank im Voraus!

[Mario A. (itse_me_mario)]

Ich hänge mal noch den Code inkl. dem Schreiben der Filterkoeffizienten 
in ein txt-File an.

[foobar (Gast)]

Hmm... scanf %d und ein int16_t?  Passt das?

Ansonsten: wenn du unerwartete Ergebnisse bekommst und schon vermutest, 
dass Overflows stattfinden, dann check das doch mal.  Ggfs den Accu 
erweitern oder saturiert arbeiten.

[foobar (Gast)]

> fscanf(fp_d,"%d\n", &s16_x[0]);
> ...
> s16_prod1 = (int16_t)((s16_x[i]  * s16_e   )>>15);

Nur eine der beiden Zeile kann korrekt sein. Wenn ints 16 Bit sind, die 
erste, wenn 32 dann die zweite.

[Mario A. (itse_me_mario)]

Das könnte ich noch mit
fscanf(fp_x,"%16d\n", &s16_x[0]);
fscanf(fp_d,"%16d\n", &s16_d);
auf 16 Stellen begrenzen. Sollte aber auch ohne gehen oder?

Hmm wenn ich

s16_prod1 = (int16_t)((s16_x[i]  * s16_e   )>>15);

in

s32_Wert1 = s16_x[i]  * s16_e;
s16_prod1 = (int16_t)(s32_Wert1>>15);

ändere? Bringt das was? Ich probiers mal. Ich bin mir nicht ganz sicher 
was du meinst. Also über die Filepointer fp_x und fp_d kann ich die 
Zahlen als int16 einlesen.

Wenn ich eine Zahl im Q1.15 Datenformat mit einer anderen Zahl im Q1.15 
Datenformat multipliziere kommt eine Zahl im Q2.30 Datenformat raus. Die 
Zahl enthält dann 2 redundante Vorzeichenbits. Durch (int16_t)(Zahl>>15) 
wird daraus wieder eine Zahl im Q1.15 Datenformat, wenn ich mich nicht 
irre.
s16_x und s16_d sollten Zahlen im Q1.15 Datenformat sein.

Im Anhang mal ein Beispiel, was ich meine und die zugehörige 
Konsolenausgabe. Edit: Gerade gesehen, dass in der letzten Zeile nicht 
int32 von 0,5*0,5*0,5 sondern nur von 0,5*0,5 ausgegeben wird mit 
printf. Der Q-Wert stimmt aber.

[foobar (Gast)]

> fscanf(fp_x,"%16d\n", &s16_x[0]);

Nein, %hd oder SCNd16:
  fscanf(fp_x, SCNd16"\n", &s16_x[0]);

Alternativ:
  int tmp;
  ...
  fscanf(fp_x,"%d\n", &tmp); s16_x[0] = (int16_t)tmp;


> s32_Wert1 = s16_x[i]  * s16_e;
> s16_prod1 = (int16_t)(s32_Wert1>>15);

Bei 16-Bit-ints ist die Multiplikation nur 16 Bit.  Sie muß 32-Bit sein, 
also mind einen Operand erweitern:
   (int32_t)s16_x[i] * s16_e

[Mario A. (itse_me_mario)]

Ahh danke für den Tipp mit %hd! Sehr guter Tipp, danke!

>>Bei 16-Bit-ints ist die Multiplikation nur 16 Bit.
Ja das hab ich vorhin übersehen! Hatte ich schon ausgebessert aber danke 
für den Tipp!

Mein Code sieht jetzt so aus:

...
#define N  32             // number of filter coeffs/weights
#define NUM_ITERS  32000  // number of iterations

int16_t s16_beta = 327;         // learning rate 0.01*2^15 = 327.68
int16_t s16_w[N] = {0};       // adaptive filter weights
int16_t s16_x[N] = {0};       // adaptive filter delay line

int32_t i = 0;
int32_t n = 0;

int16_t s16_d = 0;
int16_t s16_y = 0;
int16_t s16_e = 0;

int16_t s16_prod1 = 0;
int16_t s16_prod2 = 0;

int32_t s32_acc = 0;
int32_t s32_Wert1 = 0;
int32_t s32_Wert2 = 0;

int main()
{
    ...

    else
    {
        for(n=0; n<NUM_ITERS; n++)
        {
            fscanf(fp_x,"%hd\n", &s16_x[0]); 
            fscanf(fp_d,"%hd\n", &s16_d);

            s32_acc = 0;     

            for (i = 0; i < N; i++)
            {
                //Achtung hier könnte Overflow auftreten
                s32_acc = s32_acc + (int32_t)s16_w[i] * (int32_t)s16_x[i]; 
            }

            s16_y = (int16_t)(s32_acc >> 15);
            s16_e = s16_d - s16_y;

            for (i = N-1; i >= 0; i--)
            {
                s32_Wert1 = s16_x[i]  * s16_e;
                s16_prod1 = (int16_t)(s32_Wert1>>15);

                s32_Wert2 = s16_prod1 * s16_beta;
                s16_prod2 = (int16_t)(s32_Wert2>>15);

                s16_w[i] = s16_w[i] + s16_prod2; //Achtung hier könnte Overflow auftreten

                if(i!=0)
                {
                    s16_x[i] = s16_x[i-1];  // shift data in delay line
                }

                if(i>0)
                {
                    fprintf(fp_w, "%d\n", s16_w[i]); //Filterkoeffizienten
                }
                else
                {
                    fprintf(fp_w, "%d\n", s16_w[i]);
                    fprintf(fp_w, "%dte Koeffizienten\n",n);
                }
            }
...

Ich hab das mal geplotted (siehe Anhang). Sieht jetzt schon mal viel 
besser aus. Im Vergleich zur Double-Implementierung ist noch Luft nach 
oben. Vllt muss ich mal eine höhere Lernrate testen. Ich hab mal 
gelesen, dass ein zu kleines beta bei Festkommarealisierungen sonst dazu 
führt, dass sich die Koeffizienten ab einem gewissen Punkt nicht mehr 
ändern.

Vllt sollte ich es auch mal mit 24 Bit versuchen (also Q0.23).

Vielen Dank foobar! Der Tipp mit %hd hat mir echt geholfen! :)


Edit: Habe mal einen Vergleich von Double und Q15 angehängt (beide mit 
beta=0,2).

[Mario A. (itse_me_mario)]

Ich hätte jetzt mal versucht das ganze auf den NLMS Algorithmus 
anzuwenden, bevor ich den LMS weiter optimiere.

An sich kommt nur eine Schleife

for(j = 0; j < N; j++)
{
   s32_input_signal_power = s32_input_signal_power + ...
   (int32_t)s16_x[j] * (int32_t)s16_x[j];
}
s16_input_signal_power = (int16_t)(s32_input_signal_power>>15);

und eine zusätzliche Division durch den ermittelten Wert 
s16_input_signal_power dazu:

for(i = N-1; i >= 0; i--)
{ 
// w = w + (beta*e*x)/(x*x)

s32_Wert1 = s16_x[i]  * s16_e;
s16_prod1 = (int16_t)(s32_Wert1>>15);

s32_Wert2 = s16_prod1 * s16_beta;
s16_prod2 = (int16_t)(s32_Wert2>>15);

s16_div = s16_prod2 / s16_input_signal_power;
s16_w[i] = s16_w[i] + s16_div; //Achtung hier könnte Overflow auftreten
...
}

Der Rest bleibt gleich. So funktioniert das ganze leider noch nicht. Ich 
habe auch mal versucht, dass ich s16_prod2 vor der Division mit 
verschiedenen Werten multipliziere bzw. einen größeren Datentyp für den 
Zähler, aber das hat leider auch NOCH nichts gebracht. Ich schau mir das 
morgen nochmal an. Jetzt geht grad nichts mehr. ;)

[foobar (Gast)]

Wenn du die Implementierungen vergleichen willst, solltest du beide auch 
mit gleichen Daten füttern.  Rechne mal die reduzierten int-Daten nach 
double zurück und fütter damit die Double-Version.  Dann weisst du, was 
mit den reduzierten Daten überhaupt noch erreichbar ist und was durch 
die int-Implementation zusätzlich verloren geht.

[Mario A. (itse_me_mario)]

Ich hab jetzt mal den Zähler vorher multipliziert (<<16). Anfangs sehen 
die Fehler e ganz gut aus. Irgendwann weichen die Fehler zu sehr vom 
Ideal ab (vgl. Data_E... und Data_Y... der Gleitkommarealisierung).

else
    {
        for(n=0; n<NUM_ITERS; n++)
        {
            fscanf(fp_x,"%hd\n", &s16_x[0]); 
            fscanf(fp_d,"%hd\n", &s16_d);

            s32_acc = 0;                    // compute filter output
            s32_input_signal_power = 0;

            for(j = 0; j < N; j++)
            {
                s32_input_signal_power = s32_input_signal_power + (int32_t)s16_x[j] * (int32_t)s16_x[j];
            }

            s16_input_signal_power = (int16_t)(s32_input_signal_power>>15);

            for(i = 0; i < N; i++)
            {
                s32_acc = s32_acc + (int32_t)s16_w[i] * (int32_t)s16_x[i]; 
            }

            s16_y = (int16_t)(s32_acc >> 15);
            s16_e = s16_d - s16_y;

            for(i = N-1; i >= 0; i--)
            {
                s32_Wert1 = s16_x[i]  * s16_e;
                s16_prod1 = (int16_t)(s32_Wert1>>15);

                s32_Wert2 = s16_prod1 * s16_beta;
                s16_prod2 = (int16_t)(s32_Wert2>>15);

                s32_test1 = ((int32_t)(s16_prod2)<<15);

                s32_test2 = s32_test1 / (int32_t)s16_input_signal_power;

                s16_div = (int16_t)(s32_test2); // >> 16);

                s16_w[i] = s16_w[i] + s16_div; 
...
            }

[Mario A. (itse_me_mario)]

>>Wenn du die Implementierungen vergleichen willst, solltest du beide auch
mit gleichen Daten füttern.

Danke für den Tipp. Ich habe mal eben die verwendeten Q15 Werte für X 
und D in MATLAB als Ganzzahlen eingelesen und dann durch 2^15 geteilt. 
Dann habe ich diese wieder in ein txt-file geschrieben (als Double) und 
anschließend nochmal gelesen und in MATLAB mit den Original Double 
Werten von X und D verglichen.

Es scheint, dass das nicht das Problem ist.

[foobar (Gast)]

> Es scheint, dass das nicht das Problem ist.

Das siehst du anhand des grafischen Vergleichs?  Da sind, wenn's hoch 
kommt, gerade mal die oberen 7-8 Bit zu erkennen und natürlich stimmen 
die überein - die Unterschiede sind in den abgeschnittenen unteren 32-48 
Bits.

[Mario A. (itse_me_mario)]

Nein ich hab mir auch die Abweichung plotten lassen und die Differenz in 
den Arrays überprüft. Aber so vom Prinzip her hast du es gemeint oder?

Ich hab jetzt mal zusätzlich noch
if(s16_input_signal_power == 0){ s16_input_signal_power = 1; }
eingebaut. Bin weiter am schauen woran´s liegt.

Bei der Division geht wohl viel verloren. Vllt sollte ich für den 
Quotienten 64 Bit verwenden. Ich versuch das mal.

[Mario A. (itse_me_mario)]

Leider hab ich noch nicht die Lösung meines Problems gefunden.. Schau 
mir morgen mal die einzelnen Koeffizienten w[n] genauer an.

[Mario A. (itse_me_mario)]

Meine Koeffizienten w werden scheinbar irgendwann zu groß... Dadurch 
wird dann auch y zu groß. Ich weiß noch nicht wie ich das umgehe. Hab 
jetzt zus. noch eine Overflow/Underflowerkennung drin mit Sättigung.