FFTW - Ergebnis zu hoch



  • Hallo zusammen,

    ich wollte mich in die FFT mit fftw (www.fftw.org) einarbeiten und habe dazu ein kleines Testprogramm geschrieben.
    Zuerst erstelle ich eine Summe von 5 Sinuskurven mit der Frequenz von 10-50Hz. Amplitude ist jeweils 10.
    Nach ausgeführter FFT wird der Betrag des Ergebnisses genommen. Ich halte mich dabei an das Beispiel von http://www.fftw.org/fftw2_doc/fftw_2.html.

    double *in;
    	double *out;
    	double *power_spectrum;
    	int i, j;
    	double dt = 0.001;
    	int N = 1024;
    	double fi = 0;
    // Signal creation
    	for (i=0; i<N; i++)
    	{
    		for (j=10; j<=50; j=j+10 )
    		{
    			in[i] += 10*sin(2*M_PI*j*fi);
    		}
    		x_time[i] = fi;
    		fi += dt;
    	}
    // FFT
    	my_fft_plan = rfftw_create_plan(N, FFTW_REAL_TO_COMPLEX, FFTW_ESTIMATE);
    	rfftw_one(my_fft_plan, in, out);
    
    //Magnitude
    	power_spectrum[0] = out[0]*out[0];  /* DC component */
    	for (i=1; i<(N+1)/2; i++)
    		power_spectrum[i] = out[i]*out[i] + out[N-i]*out[N-i];
    	 if (N % 2 == 0) /* N is even */
    		power_spectrum[N/2] = out[N/2]*out[N/2];  /* Nyquist freq. */
    

    Leider habe ich keine Möglichkeit das graphische Ergebnis als Bild anzuhängen, also versuche ich es auszuformulieren:

    Es werden fünf Spektrallinien an der richtigen Stelle (in x-Richtung) angezeigt, leider sind die Amplituden viel zu hoch (um 2e+07).
    Außerdem sind die Amplituden der fünf Peaks deutlich unterschiedlich hoch (von 1e+07 bis 2.5e+07)

    Ich habe auch versucht die Beträge mit

    sqrt( sqr(RE) + sqr(IM))
    

    zu berechnen, dadurch werden die Werte zwar deutlich kleiner, liegen aber immer noch zwischen 3.000 und 5.000

    Kann mir jemand erklären was ich falsch mache?
    Vielen Dank im Voraus!!!

    Gruß zorin



  • http://www.fftw.org/fftw3_doc/The-1d-Discrete-Fourier-Transform-_0028DFT_0029.html#The-1d-Discrete-Fourier-Transform-_0028DFT_0029

    > FFTW computes an unnormalized transform, in that there is no coefficient in front of the summation in the DFT. In other words, applying the forward and then the backward transform will multiply the input by n.



  • Ok, vielen Dank für den Tip. Hab's leider trotzdem noch nicht verstanden:

    Die FFT wird ohne Faktor durchgeführt, also ohne Amplitude des Sinus-Signals?
    Wo ist dann der Zusammenhang zwischen Eingangs-(Zeit-)Amplitude und Resultat(Frequenz)-Amplitude?

    Und wie erklärt sich der Unterschied in den verschiedenen Peaks? Die sollten doch zumindest gleich (falsch) sein?!

    Gruß Zorin



  • Das könnte mit dem randeffekt zusammenhängen, weil du kein periodisches signal hast. Was passiert denn, wenn du das signal das du tranformierst länger machst?



  • Hi,

    Was meinst Du mit "kein periodisches Signal"? Das Signal besteht doch aus fünf (klinisch) reinen Sinussignalen?!

    Ok folgender Graph wurde mit folgenden Eckdaten erstellt:
    N=1024; Signalamplitude = 10;

    power_spectrum[i] = sqrt(out[i]*out[i] + out[N-i]*out[N-i]);
    

    http://s14.directupload.net/file/d/3092/yal5qjf8_png.htm

    Wenn ich N auf 4096 erhöhe kommt das bei raus:
    http://s7.directupload.net/file/d/3092/x7lees2g_png.htm

    Meinst Du das mit "länger machen"? Ich kann leider keinen Zusammenhang zwischen Eingangsamplitude, FFT-Parametern und FFT Ergebnis feststellen....

    Gruß Zorin



  • Das Problem ist, dass die FFT annimmt, dass das Signal periodisch ist. Dadurch entstehen am Rand komische Artefakte, wenn Anfang und Ende nicht 100% zusammenpassen. Stell Dir vor, Du würdest eine Kopie des Plots rechts daneben anschließen, dann entsteht dort ein kleiner Sprung, und dieser ist äußerst hochfrequent.

    Wenn Du das Signal verlängerst, dann werden diese Anteile immer hochfrequenter und fallen letztlich nicht mehr so auf. -- Das siehst Du auch in Deinem Bild, die Peaks sind im zweiten Fall deutlich schärfer und auch die Höhen liegen dichter beieinander.

    Du kannst mal versuchen die Länge so zu wählen, dass sie genau auf eine gemeinsame Periode aller Sinuskurven fällt, dann dürftest Du ein reines Spektrum bekommen. In der Praxis behilft man sich normalerweise mit einer Windowing-Funktion, die das Signal vorne und hinten differenzierbar auf 0 drückt. Dadurch entstehen zwar auch ein paar "Fehler" im Spektrum, aber die haben nicht so einen riesigen Effekt.



  • Was für eine Verteilung hast Du denn erwartet? Willst Du das Fourierspektrum? Warum berechnest Du dann die Leistungsdichte?

    Schau dir nochmal die Grundlagen an, die FFT macht im wesentlichen

    Y(k)= Sum_(i=0)^n x(i) exp(-2*pi*j/n)^((i-1)*(k-1))

    nur in schlau. Wenn Du da jetzt einen mit, sagen wir, 10000 Punkten gesampleten sin(t) t aus [0,100] reinsteckst, dann kannst Du ja kaum einen Peak mit der Höhe 1 (oder pi oder 2*pi oder 1/pi oder wie auch immer deine bevorzugte Normierung aussieht) erwarten, sondern der wird ne ganze Ecke höher sein, denn Du hast ja mehr Punkte reingenommen. Dividiere die Anzahl n mal wieder raus und dann deinen Lieblings-Skalierungsfaktor drauf und Du solltest das DFT-Spektrum haben (evtl. muß Du auch noch mit 2 oder 1/2 multiplizieren, wenn dein Frequenzspektrum nur halb- oder beidseitig sein soll; sollte hier deinen Bildern zufolge aber passen). Daraus kannst Du dann immer noch dein Leistungsdischtespektrum schätzen, wenn Du magst.

    Sieht das dann schon eher nach deiner Vorstellung aus?



  • zorin59 schrieb:

    Und wie erklärt sich der Unterschied in den verschiedenen Peaks? Die sollten doch zumindest gleich (falsch) sein?!

    Das nennt man im Fachjargon "spectral leakage". Du hast zwar ein periodisches Signal, aber die Blocklänge deiner FFT ist kein ganzzahliges Vielfaches der Periodendauer. Deswegen triffst du mit den Frequenzen auch nicht exakt die "FFT-bins" sondern bist irgendwo dazwischen. Weil du implizit ein Rechtecksfenster verwendest, wird das Spektrum mit einem Sinc verschmiert. Was du also hinterher siehst, ist eine Abtastung von überlagerten Sinc-Kurven. Und bei der Abtastung triffst du nicht extakt die Peak-Positionen.

    Wenn du an den abgetasteten Werten im Spektrum die Amplitude einer harmonischen Komponente ablesen willst und die Frequenz dieser beliebig ist, musst du das Maximum des Peaks entweder interpolieren oder du verwendest vor der FFT ein Flat-Top Fenster. Das Flat-Top-Fenster führt im Spektrum dazu, dass die Peaks flach und breit genug sind, dass man die Höhe nicht per Interpolation bestimmtn muss, sondern direkt an dem Wert ablesen kann, der am nächsten am Peak dran sitzt.

    Interpolation bekommt man übrigens fast geschenkt. Dazu kann man ein "zero padding" vor der FFT durchführen und die "FFT vergrößern". Also, statt 256 Werte in die FFT zu stecken, steckst du z.B. 512 Werte rein, wovon die letzten 256 Werte 0 sind. Die ersten 256 Werte sollten "vernünftig" "gefenstert" werden.

    Was du machen solltest, hängt wirklich davon ab, welches Problem du zu lösen versuchst. Interessieren dich die Amplituden von einzelnen harmonischen Komponenten, verwendet man typischerweise ein Flat-Top-Fenster und skaliert auch dementsprechend. Interessiert man sich für die spektrale Leistungsdichte, verwendet man andere Fenster (z.b. Blackman) und skaliert nochmal anders.

    Du musst auch berücksichtigen, dass sich ein "Sinus" aus zwei Frequenzen zusammensetzt: der positiven und der negativen. Reellwertige Signale sind ja nur Spezialfälle von komplexwertigen mit verschwindendem Imaginärteil.

    cos(f*t) = 0.5*( exp(f*i*t) + exp(-f*i*t) )      mit  -1 = i*i
                         ^            ^^
    

    Das kommt auch bei der Skalierung mit ins Spiel.



  • zorin59 schrieb:

    power_spectrum[i] = sqrt(out[i]*out[i] + out[N-i]*out[N-i]);
    

    wenn es wirklich proportional zur "power", also Leistung sein soll, dann muss die Wurzel auf der rechten Seite weg.

    Interessant ist auch noch die Frage, ob du wirklich Leistung haben willst, oder vielleicht doch die spektrale Leistungsdichte. Das skaliert man nämlich anders. Die Leistung im einem bestimmten Frequenzband ist nämlich das Integral der spektralen Leistungsdichte über die Frequenzen. Wenn du jetzt die Blöckgröße der FFT änderst, dann beeinflusst du natürlich auch die Bandbreite der schmalen Frequenzbänder, in die du das Signal per FFT zerlegst.

    Wenn du bei variabler FFT-Blockgröße und einem bestimmten Sinus-Test-Signal immer dieselbe Amplitude ausrechnen willst, dann gehst du über die Leistung.

    Wenn du bei variabler FFT-Blockgröße und einem breitbandigem Rauschen immer dieselbe Leistungsdichte ausrechnen willst, dann skalierst du entsprechend, um die Leistungsdichte zu erhalten.

    Beide Ansätze sind inkompatibel. Damit meine ich, es gibt keine einzige Skalierung, mit der du beides erhältst, also mit der du sowohl Amplituden von harmonischen Anteilen als auch die Leistungsdichte direkt ablesen kannst.

    Aber man kann relativ einfach hin und herrechnen:
    spektrale Leistungsdichte = Leistung / Bandbreite.
    wobei sich die Bandbreite aus der Abtastfrequenz und der FFT-Blockgröße ergibt. Beispiel: fs=44100Hz, N=1024, dann haste eine Bandbreite von fs/N=43Hz.



  • Ok, habe nun den Rat von Jesper und Daniel befolgt und das Abtastintervall als Vielfaches der Frequenzen gesetzt sowie den Output durch N/2 geteilt. Nun kommen saubere Peaks raus.

    Mir war nicht bewusst, dass das Signal am Ende das Abtastintervall so einen starken Einfluss hat. Muss mich also mit dem Thema "Fenster" etc. genauer auseinandersetzen. Meine Vorstellung auch unter realen Bedingungen einen sauberen Peak (mit korrekter Höhe) zu bekommen war wohl zu optimistisch.

    Eigentlich ist der Variablenname "power_Spectrum" falsch bezeichnet. Ich ging von einer Eingangsspannung aus.

    Vielen Dank für Eure Hilfe!

    Gruß Zorin




Anmelden zum Antworten