Falsche Amplituden aus DFT?

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.

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..

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.

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

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.

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!

Dann rechne doch mal c_k per Hand aus.

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?

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.

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:

1

clear all;per=2.22;dt=0.005;nstart=round(per/dt);

2

%Gerades n erzeugen

3

if mod(nstart,2)

4

    n=nstart+1;

5

else

6

    n=nstart;

7

end

8

9

%Frequenzauflösung

10

df=1/(n*dt);

11

%Samplingfrequenz fs

12

fs=1/dt;

13

%Grenze zur Beschränkung der Plots auf interessanten Bereich

14

grz=100;

15

%Diskretisierung

16

for k=1:1:n

17

    t(k)=dt*(k-1);

18

    y(k)=2*sin(2*pi*t(k))+sin(12*pi*t(k))+cos(12*pi*t(k));

19

end

20

xa=[0:fs/n:fs-1/n];

21

y=y';

22

23

%Plot der kontinuierlichen und diskretisierten Funktion

24

subplot(2,3,1), plot(t,y);grid on

25

xlabel('t[s]')

26

ylabel('f(t)')

27

title('kontinuierliches Signal');

28

29

subplot(2,3,4), plot(xa,y,'.'); grid on

30

xlabel('dt*i')

31

ylabel('y(i)')

32

title('Abtastwerte, diskret');

33

34

%DFT Algorithmus 

35

36

ft=dftmtx(n);

37

%Spektralkoeffizienten c erzeugen

38

c_dft=(1/n)*ft*y;

39

%Rundung

40

%c_dft=round(c_dft*100)/100;

41

%Betrag der Koeffizienten

42

amp_dft=zeros([1,n]); ph_dft=zeros([1,n]);

43

44

%Verdopplung der Amplitude, dadurch 

45

%Berücksichtiung der positiven und negativen Frequenzen

46

%des Spektrums

47

amp_dft=2*abs(c_dft);

48

49

%Phase DFT, Minus für Berücksichtigung unterschiedlicher Konvention der

50

%Phasen

51

ph_dft=-angle(c_dft);

52

53

%Plot auf interessanten Frequenzbereich beschränken 

54

xa(grz:end)=[];

55

amp_dft(grz:end)=[];

56

ph_dft(grz:end)=[];

57

58

%Plot der einseitigen Spektren

59

subplot(2,3,2), plot(xa,amp_dft,'-b*');grid on

60

xlabel('n');

61

ylabel('abs(c_d_f_t)');

62

title('diskr. DFT Amplitudenspektrum');

63

ylim([0,2.25]);

64

%set(gca,'YTick',-10:0.25:10)

65

%set(gca,'XTick',-100:1:100)

66

67

subplot(2,3,5), plot(xa,ph_dft,'-b*');grid on

68

xlabel('n');

69

ylabel('phi [rad]');

70

title('diskr. DFT Phasenspektrum');

71

%set(gca,'YTick',-20:0.4:20)

72

%set(gca,'XTick',-20:1:20)

73

74

75

%FFT Algorithmus 

76

77

c_fft=zeros([1,n]);av=zeros([1,0.5*n]);bv=zeros([1,0.5*n]);

78

diagmx=zeros(0.5*n,0.5*n);a_fft=zeros([1,n]);b_fft=zeros([1,n]);

79

ph_fft=zeros([1,n]);

80

81

%Aufspalten von y in gerades a, ungerades b

82

for k=1:0.5*n

83

    av(k)=y(2*k);

84

    bv(k)=y(2*k-1);

85

end

86

av=av';

87

bv=bv';

88

89

%Berechnung von g, u

90

fftmx=conj(dftmtx(0.5*n));

91

92

%Diagonalmatrix für b Multiplikation

93

for k=0:0.5*n-1

94

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

95

end

96

97

%Multiplikationen:

98

m_a=fftmx*av; 

99

m_b=diagmx*fftmx*bv;

100

101

%Obere und untere Hälfte von c erzeugen:

102

gv=(1/n)*(m_a+m_b);

103

uv=(1/n)*(m_a-m_b);

104

105

%Spektralkoeffizienten aus gv und uv zusammensetzen

106

for k=1:n

107

    if k<=0.5*n

108

        c_fft(k)=gv(k);

109

    else

110

        c_fft(k)=uv(k-0.5*n);

111

    end  

112

end

113

114

%Amplitudenspektrum aus FFT

115

amp_fft=abs(c_fft);

116

117

%Phasenspektrum aus FFT

118

%Runden

119

%c_fft=round(c_fft*100)/100;

120

ph_fft=angle(c_fft);

121

%FFT-Spektrum (einseitig) auf interessanten Bereich beschränken

122

amp_fft(grz:end)=[];

123

ph_fft(grz:end)=[];

124

125

%Plot der einseitigen Spektren

126

subplot(2,3,3),plot(xa,2*amp_fft,'-b*');grid on

127

title('diskr. FFT Amplitdenspektrum');

128

xlabel('n');

129

ylabel('abs(c_f_f_t)');

130

%set(gca,'YTick',-10:0.25:10)

131

%ylim([0,2.25]);

132

133

subplot(2,3,6),plot(xa,ph_fft,'-b*');grid on

134

title('diskr. FFT Phasenspektrum');

135

xlabel('n');

136

ylabel('phi [rad]');

137

%set(gca,'YTick',-20:0.2:20)

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

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?

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?

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.

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

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)

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!

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.

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!

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.