Fourier Matlab

[Dennis (Gast)]

Hallo,
ich habe ein Verständnis Problem mit meiner Ausgabe. Ich erhalte zwei 
Spektren, wenn ich eine Sinus Funktion mit der FFT transformiere.
Hier mein Programm.
Vielleicht kann mir jemand erklären, wieso ich ein weiteres Spektrum bei 
900Hz erhalte. (Bzw. bei Abtast-f)
Hier mein Programm:

clc

%%%%%%%%%%%%%%%%%%%
%  Eingabe:       %
%%%%%%%%%%%%%%%%%%%
a = 1;
w = 2*pi;
f = 100;
n = 0.001;
time = (0:n:1)';
time1 = (n:n:1)';
abtast = floor(1/n);
lt = length(time);


%%%%%%%%%%%%%%%%%%%
%  Schwingung     %
%%%%%%%%%%%%%%%%%%%
v_sin = a.*sin(w.*f.*time);
subplot(2,1,1)
plot(time, v_sin)
title('Schwingung')
xlabel('Zeit [s]')
ylabel('Amplitude')
grid on;


%%%%%%%%%%%%%%%%%%%
%  Fourier        %
%%%%%%%%%%%%%%%%%%%
z = fft(v_sin,abtast)
Leistungsspektrum = z.*conj(z)/abtast; %Leistungsdichtespektrum
b = length(Leistungsspektrum);
x = (0:((abtast-1)));
subplot(2,1,2)
stem(x,Leistungsspektrum(1:abtast));
title('Spektrum')
xlabel('Frequenz [Hz]')
ylabel('Amplitude')
grid on;

%%%%%%%%%%%%%%%%%%%%
%  Ausgabe         %
%%%%%%%%%%%%%%%%%%%%
[MaxLeistungsspektrum,freq] = max(Leistungsspektrum);
freqMax = (freq-1);
MaxAmplitude = MaxLeistungsspektrum
MaxFrquenz = freqMax

[die ??? (Gast)]

Dennis wrote:
> Ich erhalte zwei Spektren, wenn ich eine Sinus Funktion mit der FFT
> transformiere.

Nicht ganz. Du bekommst ein Spektrum geliefert. Das Spektrum spiegelt 
sich bei der halben Abtastfrequenz (Nyquistfrequenz).

  http://de.wikipedia.org/wiki/Nyquist-Frequenz


> Vielleicht kann mir jemand erklären, wieso ich ein weiteres Spektrum bei
> 900Hz erhalte.

Vorsicht, nicht die Begriffe durcheinander schmeißen. Du meinst 
sicherlich nicht Spektrum sondern Spektrallinie.

[Beate (Gast)]

Hallo,

ich habe ein ähnliches Programm wie der Dennis geschrieben.
Funktioniert auch alles super.
Ich habe allerdings wohl eher ein mathematisches Verständnis Problem.
Vielleicht kann mir trotzdem jemand behilflich sein.
Ich erzeuge einen Sinus mit Amplitude=1 und einer Frequenz von 10Hz
Ich taste nun das Signal mit 1000 ab. Erhalte einen schönen Peak.
Mein Problem ist nun die Amplitude des Spektrums.
In meinem Fall erhalte ich eine Amplitude von genau 250.
Erhöhe ich die Ampltude des Sinus auf 2 erhöht sich die Amplitude des 
Spektrums auf 1000.
Wie kommen die Amplituden im Spektrum zustande? Gibt es da einen 
Mathematischen Hintergrund?

Ich hoffe mir kann jemand helfen.
Danke und Liebe Grüße
Beate

[gast (Gast)]

da im spektrum die leistung angegeben wird. doppelte amplitude -> 
vierfache leistung

[Beate (Gast)]

Ok, aber so ganz habe ich es leider noch nicht verstanden.
Ich dachte immer, dass ein Spektrum eine Abbildung der Frequenz 
(x-Achse) x Amplitude (y-Achse) von Sinusoiden ist.

So hätte ich ja dann eine Amplitude von 1 erhalten müssen und nicht von 
250. Gibt es denn eine Formel für die Leistungsberechnung.

Danke.
LG Beate

[Patric (Gast)]

Bevor ihr hier mit Leistungsdichtespektrum usw. anfangt, versucht 
erstmal die DFT(FFT) zu verstehen. Mit fft(x) ist es leider nicht getan. 
Ich hab euch mal ein Matlab-Programm zur Demonstration angehängt. Es 
sind einige Signalbeispiele enthalten, die bestimmte Eigenschaften der 
DFT darstellen. Versucht das mal zu verstehen und meldet euch bei Bedarf 
hier wieder.

%
% Berechnung der FFT und Darstellung in Bins und Frequenzen
%
clear all;

% Parameter
N = 1000;
fs = 1e3;
f = (0:(N-1))./N*fs;

% numerische Fehler vermeiden
err = 1e-6;

% Zeit
t = 0:1/fs:(N-1)*1/fs;

% Signal erzeugen
x = 5 + cos(2*pi*200.*t);
%x = 5 + cos(2*pi*1000.*t);
%x = -5 + cos(2*pi*1000.*t);
%x = sin(2*pi*500.*t+pi/2);
%x = sin(2*pi*500.*t+pi/4);
%x = 10*cos(2*pi*1000.*t) + 8*cos(2*pi*200.5.*t);
%x = cos(2*pi*2000.*t);

% FFT
X = fft(x,N);

% Betrag
betrag = abs(X)/N;
% Winkel
winkel = angle(X);
for p = 1:N
    if (abs(real(X(p))) < err & abs(imag(X(p))) < err)
        winkel(p) = 0;
    end;
end;


figure(1);
subplot(221), stem(0:(N-1), betrag); grid on;
title('abs X(f)'); xlabel('Bins');
subplot(223), stem(f, betrag); grid on;
title('abs X(f)'); xlabel('Frequenz');

subplot(222), stem(0:(N-1), winkel); grid on;
title('angle X(f)'); xlabel('Bins');
subplot(224), stem(f, winkel); grid on;
title('angle X(f)'); xlabel('Frequenz');

[Sebastian (Gast)]

Wie Patrick in seinem Programm schon richtig zeigte, muss die Skalierung 
bei fast allen Signalen mit dem Faktor 1/N stattfinden.

betrag = abs(X)/N;

Allerdings ergibt dann der Betrag eines Sinussignals der Amplitude 1 
(0.5 PP) im Frequenzspektrum 0.5 bzw. -6db. Genau diese Tatsache habe 
ich bisher aber nicht verstanden...wieso die Halbierung?

Außerdem liefert die Skalierung mit 1/N bei einem Diracimpuls ein 
falsches Ergebnis. Auch hier konnte tappe ich bisher im Dunkeln...

[joep (Gast)]

Du hast ja zwei Spektrallinien (einmal im positiven und einmal im 
negativen Frequenzbereich) mit jeweils Amplitude von 0.5 => 2*0.5=1.

Was liefert die Skalierung mit 1/N bei einem Dirac für ein Ergebnis und 
was für eins erwartest du?

[Patric (Gast)]

> Allerdings ergibt dann der Betrag eines Sinussignals der Amplitude 1
> (0.5 PP) im Frequenzspektrum 0.5 bzw. -6db. Genau diese Tatsache habe
> ich bisher aber nicht verstanden...wieso die Halbierung?

Du solltest mein Programm nochmal etwas genauer anschauen. Die 
Beispielsignale mit Phase sind hier sehr interessant.

Die DFT liefert immer zwei Ergebnisse. Eines bei der Frequenz f und ein 
komplex konjugiertes bei fs-f. Beide ergeben vektoriell addiert wieder 
den Betrag (Amplitude) des Zeitsignals. D.h. die -6dB, die du 
beobachtest kommen nur dann zustande wenn die Phase des Zeitsignals so 
liegt, dass die Addition der Amplituden mit der vektoriellen Addition 
übereinstimmt. Das sieht man in dem Beispielprogramm sehr deutlich.

Thema Dirac:
Das Programm oben liefert genau das richtige Ergebnis. Wenn man z.B.

% Dirac erzeugen
x = 0.*t; x(floor((N-1)/2)) = 1;

ergänzt. Erhält man im Frequenzbereich eine Konstante bei 10^-3. Alle 
Komponenten der FFT addiert (N=1000) ergibt 1.

[Sebastian (Gast)]

Auch wenn es schon einige Zeit her ist, möchte ich dir für deine 
schnelle Antwort danken...es hat mir sehr weitergeholfen. Ich habe 
selbst in C einen FFT Algorithmus erstellt und nicht deine Matlabfile 
benutzt. Ansonsten hätte ich das wohl eher nachvollziehen können.

Nun habe ich allerdings ein neues Problem...ich befasse mich zwar nicht 
mit Matlab aber immerhin geht es um die Fourieranalyse.
Ich habe einen FFT, DFT, IFFT und IDFT Algorithmus programmiert doch 
leider liefert die inverse DFT kein korrektes Ergebnis.

Laut Formel der IDFT muss man zunächst den konjugiert komplexen Betrag 
aus den Spektrallinien der DFT bilden. Nun kann man doch denselben 
Algorithmus der DFT für die Rücktransformation verwenden...abgesehen vom 
Faktor 1/N der noch hinzukommt?? Leider entspricht das Zeitsignal aber 
in keinster Weise mehr dem ursprunglichen Signal vor der Transformation 
in den
Frequenzbereich.

Entscheident ist wohl auch noch, ob es sich um ein reelles 
Signal...falls nicht, müßte man nach der IDFT nochmal konjugiert 
komplexe Werte bilden. Oder habe ich da wieder etwas falsch verstanden.

Gruß Sebastian

[T. H. (pumpkin)]

Sebastian wrote:
> Laut Formel der IDFT muss man zunächst den konjugiert komplexen Betrag
> aus den Spektrallinien der DFT bilden.

Sicher?

> Nun kann man doch denselben
> Algorithmus der DFT für die Rücktransformation verwenden...abgesehen vom
> Faktor 1/N der noch hinzukommt??

Und dem nicht negativen Exponenten. Macht aber nix. (Ob sich der Zeiger 
links oder rechts rum dreht ist egal, man hat nur die Konjugierte - was 
man aber bedenken muss.)

> Entscheident ist wohl auch noch, ob es sich um ein reelles
> Signal...falls nicht, müßte man nach der IDFT nochmal konjugiert
> komplexe Werte bilden. Oder habe ich da wieder etwas falsch verstanden.

Nein, wieso? Die DFT verhält sich "gleich" für reelle und komplexe 
Signale.

[Sebastian (Gast)]

Den konjugiert komplexen Betrag muss ich bilden, um den Algorithmus der 
DFT für die Rücktransformation verwenden zu können.
Ansonsten müsste man statt...

cos(x)-j*sin(x) nämlich cos(x)+j*sin(x) rechnen. Das ist der einzige 
Grund weshalb der Imaginärteil vor der IDFT mit -1 multipliziert wird.

Ich habe nun aber die Lösung gefunden, aber ich verstehe sie nicht.
Nach dem die IDFT Real und Imaginärteil berechnet hat, muss ich sie von 
einander abziehen, damit das richtige Zeitsignal herauskommt.

Aber weshalb die Subtraktion???

[Günter -. (guenter)]

Sebastian wrote:
> Den konjugiert komplexen Betrag muss ich bilden, um den Algorithmus der
> DFT für die Rücktransformation verwenden zu können.
> Ansonsten müsste man statt...

Meines Wissens ist die Berechnung:

IFFT = Conj( FFT( Conj(<input data>) )

[Sebastian (Gast)]

Da hast du auch recht...

Aber bei einem reellen Zeitsignal das zunächst in den Frequenzbereich 
transformiert wurde, kann nach der Rücktransformation der letzte Schritt 
auch weggelassen werden (also die erneute Bildung des konjugiert 
komplexen Betrages). Würde man ihn aber dennoch ausführen, müßte man 
dann ja Real- und Imaginärteil addieren.

Es das überhaupt der richtige Weg um das Zeitsignal zu erhalten..also 
die Addition von RT und IT??