Falsche Amplituden aus DFT?

[Janek U. (patalmqx)]

Hallo Freunde der Mathematik,

ich beschäftige mich gerade mich der Fourier-Transformation. Weiter 
möchte ich ein Beispiel für die DFT und FFT durchrechnen, aber bei der 
DFT scheint ein Faktor bei meinem Amplitudenspektrum gegenüber dem aus 
der analytischen Reihenentwicklung nicht zu stimmen.
Ich habe die Frage auch im gomatlab Forum gestellt. Dort sind auch die 
weiteren Informationen inkl. dem matlab-DFT code zu finden:

http://www.gomatlab.de/falsche-amplitude-dft-t28657.html

Habt Ihr eine Idee, wieso die Amplituden der DFT nur halb so groß sind? 
Mit welcher Begründung/Erklärung könnte ich dort noch einen Faktor 
einführen?

Vielen Dank schonmal!

[Dumdi D. (dumdidum)]

Deine Rechnung ist falsch. Das erwartete Ergebnis ist zu hoch, da Du die 
'negativen Frequenzen' wegschmeisst. Genauer:

cos(x)=1/2 (exp(ix)+exp(-x))

d.h. bei der positiven Frequenz wuerdest Du bei der DFT von cos dann nur 
1/2 erwarten.

[Janek U. (patalmqx)]

Danke für die Antwort(beim zweiten exp fehlt noch ein "i" nehme ich 
an?)! Also die Handrechnung, die zu den Amplituden 2 und Sqrt(2) geführt 
hat, ist falsch? Und das DFT Ergebnis demnach richtig?
Nicht, dass ich den Fehler in der DFT suche und die manuelle 
Fourierkoeffizientenberechnung der wahre Übeltäter ist..

[Dumdi D. (dumdidum)]

mathias ulf schrieb:
> beim zweiten exp fehlt noch ein "i" nehme ich an?)
ja

> Also die Handrechnung, die zu den Amplituden 2 und Sqrt(2) geführt
> hat, ist falsch?
auch ja

> Und das DFT Ergebnis demnach richtig?

und nochmals ja.

[Janek U. (patalmqx)]

Danke nochmals für die Antwort. Ich habe meine Handrechnung nochmal 
durchgespielt und kann leider keinen Fehler entdecken. Meine falschen 
(?) Koeffizienten b1,b3,a3 würden beim Einsetzen in die Fourierreihe 
anscheinend zur gewünschten Ausgangsfunktion führen. Die Formel 
An=Sqrt(an^2+bn^2) ist auch richtig?

Hier die Rechnung:
http://www.imgbox.de/show/img/svekMQlTNk.PNG
http://www.imgbox.de/show/img/OifsrxAa2e.PNG

[Dumdi D. (dumdidum)]

Ah, endlich schreibst Du mal wie Du gerechnet hast. Kurze Antwort : Das 
sind zwar Fourier-Integrale, aber keine DFT was Du da machst. Siehe auch
http://en.wikipedia.org/wiki/Discrete_Fourier_transform

gut, Du kannst jetzt sagen, eine DFT sollte in zumindest eine Naeherung 
zu den Fourierintegralen sein, also sollte auch der selbe Vorfaktor 
herauskommen. Dazu : Was Deine Funktion (richtig) ausrechnet ist die 
komplexe Fourierreihe: http://de.wikipedia.org/wiki/Fourierreihe
Den Zusammenhang zu den reellen Komponenten steht auc bei Wiki.

[Janek U. (patalmqx)]

Mit der Handrechnung wollte ich tatsächlich keine DFT durchführen. Ich 
wollte nur das Amplitudenspektrum von f(t) über die Koeff. der 
Fourier-Reihe gewinnen.
Mein Weg zur DFT (nach dem Vachenauer) ist angehangen. Dort wird 
natürlich eine komplexe Reihe genutzt.

In der englischen Wiki heißt es ja unter "Definition", dass die 
Amplitude, je nach Quelle, manchmal normiert wird (im Wiki für die DFT 
anscheinend mit dem Faktor 1). Ich habe in Matlab ja auch nur den Betrag 
des komplexen Koeffizienten gebildet, was dem entsprechen müsste. Kann 
es sein, dass mit den DFT Formeln aus dem Vachenauer eine andere 
Normierung nötig ist?

Danke erstmal für die Mühle, die Ihr euch macht!

[Dumdi D. (dumdidum)]

Dann rechne doch mal c_k per Hand aus.

[Janek U. (patalmqx)]

Ich glaube, ich weiß wo mein "Fehler" lag. Die komplexe Darstellung 
erzeugt anscheinend für reale Eingangsgrößen ein zweiseitiges Spektrum. 
Mein DFT code hat ja auch zu Beginn und Ende des Amplituden-arrays die 
Amplituden zu den positiven und negativen Frequenzen erzeugt. Also 
müsste ich die Beträge/Amplituden der positiven und negativen 
Frequenzbereiche berücksichtigen (addieren?). Dies führt dann auf die 
gewünschte Verdopplung der Amplitude.

Stimmts?

[Forist (Gast)]

mathias ulf schrieb:
> Hier die Rechnung:
Und auch einmal fürs Forum ...

[Dumdi D. (dumdidum)]

mathias ulf schrieb:
> Stimmts?

es ist fuer das Verstaendnis wirklich besser es zu rechnen. Aber 
ansonsten : Im Prinzip stimmts, obwohl Du fuer die sin Terme nicht 
addieren sondern subtrahieren musst.

[Janek U. (patalmqx)]

Hallo nochmal,

ich habe mal versucht, die DFT und FFT (die einfache Vorstufe noch ohne 
Verschachtelung) in Matlab zu programmieren. Soweit sieht das auch gut 
aus, solange ich das Fenster für die Diskretisierung als ganzzahliges 
Vielfaches der Periode meiner Funktion benutze. Dann tritt kein 
Leck-Effekt auf und die Ergebnisse stimmen sehr genau mit denen aus der 
Reihenentwicklung überein.

Jetzt habe ich den Abtastzeitraum etwas verändert z.B. 1->2.22 und 
erhalte, neben dem bekannten Verschmieren der Amplitude, merkwürdige 
Phasenkurven. Ist das auch der Leckeffekt?


Edit: Ich habs auch mal mit Runden probiert, weil die Phasen zu kleinen 
Amplituden sowieso keine große Bedeutung haben. In dem Fall hilft das 
jedoch nicht (deswegen ist der Teil im Code auskommentiert).
Ich habe leider kaum Informationen darüber gefunden, was der Leck-Effekt 
mit der Phase anstellt.


Hier mal der Code, an dem ich gerade herumdoktore:

clear all;per=2.22;dt=0.005;nstart=round(per/dt);
%Gerades n erzeugen
if mod(nstart,2)
    n=nstart+1;
else
    n=nstart;
end

%Frequenzauflösung
df=1/(n*dt);
%Samplingfrequenz fs
fs=1/dt;
%Grenze zur Beschränkung der Plots auf interessanten Bereich
grz=100;
%Diskretisierung
for k=1:1:n
    t(k)=dt*(k-1);
    y(k)=2*sin(2*pi*t(k))+sin(12*pi*t(k))+cos(12*pi*t(k));
end
xa=[0:fs/n:fs-1/n];
y=y';

%Plot der kontinuierlichen und diskretisierten Funktion
subplot(2,3,1), plot(t,y);grid on
xlabel('t[s]')
ylabel('f(t)')
title('kontinuierliches Signal');

subplot(2,3,4), plot(xa,y,'.'); grid on
xlabel('dt*i')
ylabel('y(i)')
title('Abtastwerte, diskret');

%DFT Algorithmus 

ft=dftmtx(n);
%Spektralkoeffizienten c erzeugen
c_dft=(1/n)*ft*y;
%Rundung
%c_dft=round(c_dft*100)/100;
%Betrag der Koeffizienten
amp_dft=zeros([1,n]); ph_dft=zeros([1,n]);

%Verdopplung der Amplitude, dadurch 
%Berücksichtiung der positiven und negativen Frequenzen
%des Spektrums
amp_dft=2*abs(c_dft);

%Phase DFT, Minus für Berücksichtigung unterschiedlicher Konvention der
%Phasen
ph_dft=-angle(c_dft);

%Plot auf interessanten Frequenzbereich beschränken 
xa(grz:end)=[];
amp_dft(grz:end)=[];
ph_dft(grz:end)=[];

%Plot der einseitigen Spektren
subplot(2,3,2), plot(xa,amp_dft,'-b*');grid on
xlabel('n');
ylabel('abs(c_d_f_t)');
title('diskr. DFT Amplitudenspektrum');
ylim([0,2.25]);
%set(gca,'YTick',-10:0.25:10)
%set(gca,'XTick',-100:1:100)

subplot(2,3,5), plot(xa,ph_dft,'-b*');grid on
xlabel('n');
ylabel('phi [rad]');
title('diskr. DFT Phasenspektrum');
%set(gca,'YTick',-20:0.4:20)
%set(gca,'XTick',-20:1:20)


%FFT Algorithmus 

c_fft=zeros([1,n]);av=zeros([1,0.5*n]);bv=zeros([1,0.5*n]);
diagmx=zeros(0.5*n,0.5*n);a_fft=zeros([1,n]);b_fft=zeros([1,n]);
ph_fft=zeros([1,n]);

%Aufspalten von y in gerades a, ungerades b
for k=1:0.5*n
    av(k)=y(2*k);
    bv(k)=y(2*k-1);
end
av=av';
bv=bv';

%Berechnung von g, u
fftmx=conj(dftmtx(0.5*n));

%Diagonalmatrix für b Multiplikation
for k=0:0.5*n-1
  diagmx(k+1,k+1)=exp(-2*pi*i*k/(0.5*n));
end

%Multiplikationen:
m_a=fftmx*av; 
m_b=diagmx*fftmx*bv;

%Obere und untere Hälfte von c erzeugen:
gv=(1/n)*(m_a+m_b);
uv=(1/n)*(m_a-m_b);

%Spektralkoeffizienten aus gv und uv zusammensetzen
for k=1:n
    if k<=0.5*n
        c_fft(k)=gv(k);
    else
        c_fft(k)=uv(k-0.5*n);
    end  
end

%Amplitudenspektrum aus FFT
amp_fft=abs(c_fft);

%Phasenspektrum aus FFT
%Runden
%c_fft=round(c_fft*100)/100;
ph_fft=angle(c_fft);
%FFT-Spektrum (einseitig) auf interessanten Bereich beschränken
amp_fft(grz:end)=[];
ph_fft(grz:end)=[];

%Plot der einseitigen Spektren
subplot(2,3,3),plot(xa,2*amp_fft,'-b*');grid on
title('diskr. FFT Amplitdenspektrum');
xlabel('n');
ylabel('abs(c_f_f_t)');
%set(gca,'YTick',-10:0.25:10)
%ylim([0,2.25]);

subplot(2,3,6),plot(xa,ph_fft,'-b*');grid on
title('diskr. FFT Phasenspektrum');
xlabel('n');
ylabel('phi [rad]');
%set(gca,'YTick',-20:0.2:20)

[Dumdi D. (dumdidum)]

Hallo,

bei DFT und FFT muss dasselbe herauskommen. Dein Matlab Programm 
verwendet zwar das Wort 'fft' sehr oft, aber eine echte FFT sehe ich 
erstmal nicht. Die Funktion heißt 'fft'

http://www.mathworks.de/de/help/matlab/ref/fft.html

[Janek U. (patalmqx)]

Stimmt, die richtige FFT sieht ja auch eine Verschachtelung vor (der 
Grund, warum N=2^m Abtastwerte benötigt werden). Die habe ich aus 
Zeitgründen erstmal rausgelassen und nur die Symmetrie der 2m-ten 
Einheitswurzeln ausgenutzt. Meine "FFT für arme" spart dadurch nur 
relativ wenige Operationen im Vergleich zur DFT ein, sollte aber 
trotzdem laufen. Habe mich dabei an dem Vachenauer HöMa 2 orientiert.

Gibt es gute Quellen (möglichst frei zugänglich), wo man sich über den 
Einfluss des Leck-Effekts auf das Phasenspektrum schlau machen kann?

[Dumdi D. (dumdidum)]

mathias ulf schrieb:
> Leck-Effekts auf das Phasenspektrum

Es ist kein Leck-Effekt. Bei der DFT und FFT muss auch bei der Phase 
dasselbe herauskommen, ansonsten hast Du Dich verrechnet / falsch 
programmiert. Was macht denn die Matlab FFT Funktion?

[Janek U. (patalmqx)]

Hab noch etwas herumprobiert. Der Fehler scheint irgendwo in meiner FFT 
ohne Verschachtelung zu liegen. Nach einer Weite weicht dessen Phase 
doch erheblich von den beiden anderen Graphen ab. Meine DFT scheint nach 
meinen Experimenten recht gut den Ergebnissen der Matlab-Funktion zu 
folgen.
Habe c=fft(y) mit y als Vektor mit diskreten Funktionswerten genutzt.

[Janek U. (patalmqx)]

Die Skalierung des Funktionsplots (oberste Reihe im linken Bild) ist nur 
auf dem Plot so daneben..

[Dumdi D. (dumdidum)]

mathias ulf schrieb:
> Habe mich dabei an dem Vachenauer HöMa 2 orientiert

Dein Fehler ist gluecklicheweise nicht Hoe. Vergleiche einfach mal 
Deinen Code mit der Formel, und der Fehler ist sehr schnell klar (ich 
verrate ihn jetzt aber nicht)

[Janek U. (patalmqx)]

Bekomme ich noch einen Tip, wo der Knackpunkt liegt? Ich hatte versucht, 
das angehangene umzusetzen. Ist da vielleicht schon der Fehler 
enthalten?
Das Buch habe ich leider gerade nicht mehr greifbar.
Danke schonmal!

[Dumdi D. (dumdidum)]

Schreib doch einfach mal die Formeln die Du in Deinem Skript verwendet 
hast auf, ohne in Bücher zu sehen, versuch es wieder zu einer DFT 
umzuformen, und dann sollte es klar sein. Selber rechnen ist heutzutage 
nicht mehr so angesagt, oder?

Falls dann noch Schwierigkeiten bestehen, poste doch mal Deinen 
Rechenweg, und dann wird Dir auch geholfen.

[Janek U. (patalmqx)]

Ich glaube, ich habe die Fehler gefunden.
Das "conj" in fftmx=conj(dftmtx(0.5*n)); war zuviel und die alte 
Indexgeschichte (der Vektor beginnt mit dem geraden Index 0, Matlab 
beginnt mit 1) hat meine gerade/ungerade Sortierung für die Vektoren a 
und b durcheinandergeworfen. Jetzt scheint alles zu passen. Danke für 
die geduldige Hilfe!

[Dumdi D. (dumdidum)]

mathias ulf schrieb:
> Danke für die geduldige Hilfe!

gern geschehen.

Es ist immer sehr schwer fremden Code zu debuggen, ich haette ja noch 
gesagt, dass in

diagmx(k+1,k+1)=exp(-2*pi*i*k/(0.5*n));

die 0.5 weg muessen. Da es jetzt wohl funktioniert, wuerde ich mich wohl 
irren.

Dir noch viel Erfolg.