Problem mit IFFT in Matlab

[Tom (Gast)]

Hallo zusammen,

ich habe wahrscheinlich ein kleines Verständnisproblem in Matlab. Ich 
möchte eine Übertragungsfunktion (H(w) bzw. H(f)) messen und daraus 
später die Impuls- (h(t)) und Sprungantwort ermitteln. Um das vorgehen 
zu testen habe ich mir in Matlab ein Testspektrum zusammengebastelt und 
wollte dieses zurück in den Zeitbereich transformieren.

Ich habe einen Vektor X mit 1000 Werten erstellt. Alle Werte sind 0 bis 
auf einen Wert bei 150, der 1 ist. Das soll einen Dirac bei 150Hz 
darstellen, der bei einem 150Hz Sinus auf einem Spektrumanalysator zu 
sehen wäre. Wenn mit Matlab ein Zeitsignal in den Frequenzbereich 
transformiert wird, wird das Spektrum bei der halben Nyquistfrequenz 
gespiegelt. Daher habe ich meinen Vektor X so bearbeitet, dass von 1-500 
die Werte von X übernommen werden und von 500-1000 die gespiegelten 
Werte von X(1:500).

So hätte ich theoretisch das Spektrum erstellt, dass Matlab mir mit 
einer FFT von einem Sinus anzeigen würde. Wenn ich jetzt allerdings die 
IFFT durchführe, bekomme ich zwar annähernd etwas periodisches, jedoch 
ist die Frequenz nicht bei 150Hz sonder bei ca. 333Hz.

Ist meine Berechnung falsch oder sind die Achsen nicht richtig? Wäre 
nett wenn mir jemand helfen könnte. Anbei der Codeschnippsel:

fa=1000;                             % Abtastfrequenz
n=1000;                              % Samples
t=(0:(n-1))/fa;                      % Zeitachse
f=(0:(n-1))*fa/n;                    % Frequenzachse
X=zeros(1,1000);                     % Vektor mit 0en
X(150)=1;                            % Dirac bei 150
X=[X(1:500) fliplr(X(1:500))];       % bei 500Hz spiegeln
plot(f,X);                           % Spektrum kontrollieren
x=ifft(X);                           % inverse FFT
plot(t,abs(x))                       % Signal im Zeitbereich

Vielen Dank

Gruß
Tom

[Detlef _. (detlef_a)]

Die Zeile muß so lauten:
X=[X(1:501) conj(fliplr(X(2:500)))];

Das x muß rein reel werden. Vielleicht nochmal Informationen zur FFT 
sichten!

Cheers
Detlef

[Tom (Gast)]

Hallo,

danke für den Tipp!

Das hat mich weiter gebracht. Mein Spektrum, dass ich mir zur Kontrolle 
erstelle ist falsch. Ich betrachte bei der FFT den Betrag und genau 
diesen habe ich mir zusammengebastelt. Mein Spektrum ist so allerdings 
nur rein reell, ich habe die Informationen des Real- und Imaginärteils 
vernachlässigt, da ich nur immer den Betrag angeschaut habe. Das 
richtige Spektrum (für einen Sinus) hätte einen negativen imaginären 
Dirac bei -150Hz und einen positiven imaginären Dirac bei +150Hz. Wenn 
ich mir dieses Spektrum erstelle, so das die Mitte bei 500Hz liegt 
funktioniert die IFFT einigermaßen genau. Dann bekomme ich Frequenzwerte 
zwischen 140Hz und 160Hz im Zeitbereich. Ich denke das sollte so 
stimmen, oder?!

Ich habe übrigens einen kleinen Fehler entdeckt, ich habe oben gesagt, 
dass bei der halben Nyquistfrequenz gespiegelt wird. Diese Aussage ist 
falsch. Es wird bei der Nyquistfrequenz bzw. bei der halben 
Abtastfrequenz gespiegelt.

Danke für die Hilfe

Gruß
Tom

[Detlef _. (detlef_a)]

Bitte, bitte.

>>die IFFT einigermaßen genau

Das Ergebnis der IFFT sollte rein reel sein, bis auf die Rundungsfehler, 
Größenordnung so 10^-15. Du kannst aber zwecks Erkenntnisgewinn 
umgekehrt vorgehen: Bastel einen 150Hz Sinus im Zeitbereich und kuck Dir 
dessen FFT an.

Cheers
Detlef

[Alexander L. (lippi2000)]

Man beachte immer:

reeles Signal    --> komplexes Spektrum

reeles Spektrum  --> komplexes Signal

P.S. Dein Spektrum wird mit "fftshift" zentriert, dadurch erhälst du die 
Darstellung  -N/2 < n < +N/2

[hierunddort (Gast)]

Hallo, ich hab hierzu eine Frage:

Klar ist, dass:
Ein Sinus-Signal im Zeitbereich mit der Amplitude 1 z.B. 
(x=sin(2*pi*f_signal*t)) ergibt nach FFT(x,N) ein Ergebnis im 
Zweiseitenband mit halber Amplituden. wird dieses mit
X=[X(1) 2*X(2:N/2)]           %N= Länge der FFT
auf das Einseitenband reduziert (mit der Länge der FFT normiert) 
entspricht die Amplitude des Spektrums der des Signals im Zeitbereich.

Fall 1:
Im Frequenzbereich ist ein ein Signal Y mit einem Dirac bei 150 gegeben.
Davon wird mit ifft(Y) direkt die iFFT berechnet.
Zur Amplituden normierung wird mit der Länge der iFFT multipliziert.
Dies ergibt im Zeitbereich ein Signal y(t) mit der Amplitude 1.

Wird das Signal im Spektralbereich mit
Y'=[Y(1:N/2+1) conj(fliplr(Y(2:N/2)))] umgerechnet,
ergibt dies letztendlich ein Signal y(t) mit der Amplitude 2.

Müsste zur Beibehaltung der Energie des Signals die Erweiterung des
Signals nicht etwa
Y''=[Y(1) 0.5*Y(2:N/2+1) 0.5*conj(fliplr(Y(2:N/2)))] sein?
Dies entspricht auch der oben genannten Formel für X.

Dann wäre letztendlich das Ergebnis von ifft(Y)=ifft(Y'')=|=ifft(Y')

Fall 2:
Wird allerdings im Frequenzbereich ein Signal Z mit dem Wert 1 bei allen 
Frequenzen eingesetzt Z(f)=1, dessen iFFT theoretisch einen Dirac
ergibt, so ergibt sich mit den Formeln
Z'=[Z(1:N/2+1) conj(fliplr(Z(2:N/2)))]
bzw
Z''=[Z(1) 0.5*Z(2:N/2+1) 0.5*conj(fliplr(Z(2:N/2)))]
dass Z=Z'=|=Z'' und im Zeitbereich
z(t)=z'(t)=1 aber
z''(t)=0.5  ist um den Faktor 2 zu klein.

wird mit ifft(Z) direkt die iFFT berechnet, so stimmt die Amplitude 
wieder.

Fazit:
Ich verstehe es nicht,
und wenn ich direkt die iFFT anwende, egal ob Einseitenband oder 
Zweiseitenband, komme ich zu den von mir erwarteten Ergebnissen.

[Alexander L. (lippi2000)]

Du hast schon vollkommen Recht. Die Amplitude muss ja auch doppelt so 
groß sein.

Schau dir mal die normale Fourier-Transformation an (also nicht 
DFT/FFT).

Ein Sinussignal der Amplitude 1 und Frequenz f0 ergibt eine 
Spektrallinie bei +f0 mit der Amplitude PI und eine Spektrallinie bei 
-f0 mit der Amplitude -j*pi.

Da der Umrechnungsfaktor zwischen Zeitbereich und Frequenzbereich aber 
2*pi teilt sich ebend die Energie auf beide auf.

Also eigentlich müsste die Amplitude bei +f0 mit 2*pi/2 angegeben werden 
:-)

[H. D. (hierunddort)]

Gut,
dann sind wir uns ja einig über die Eigenschaften bei einzelnen 
Sinus-Signalen, bzw. einzelnen Frequenzen.
Es wird davon ausgegangen das die Frequenz des Signals < Fs/2 ist., dann 
wird entsprechend von f(0:N/2-1) auf f(N/2:N) erweitert.

Wie sieht das aber bei Dirac-Impulsen im Zeitberiech bzw. breitband 
Signalen im Frequenzbereich aus?

Wenn das Signal meinetwegen die Amplitude 2*pi/2 von f(0) bis f(Fs) hat?
Wird dann, wie oben, das Signal von f(0:N/2-1) genommen und auf f(N/2:N) 
erweitert?
Dann gibt es den Amplitudenfehler von 2.

[Alexander L. (lippi2000)]

Nein! Das ist genau das Gleiche.

Reelle Signale haben immer ein symmetrisches Spektrum, da kommst du 
nicht herum!

Ein Dirac im Zeitbereich besitzt ein gleichbleibendes Spektrum im 
Bereich von

$$-f( \infty )\dots +f(\infty )$$

Dein Signal hat dann IMMER ein Spektrum von -fs/2 bis +fs/2.

Willst du dir demnach ein Spektrum vorgeben, dann muss jede 
Spektralseite nur die halbe Amplitude besitzen.

Bitte überprüfe den folgenden Matlab - Code:

N = zeros(1,1000);     
%Kosinusignal mit 10 Perioden in 1000 Samples 
%Achtung N(0) ist der DC-Wert
N(11) = 1;
%Jetzt muss der symmertrische Teil erzeugt werden
N = [N(1:501) conj(fliplr(N(2:500)))];
subplot(2,1,1);
plot(N);          %Hier geht auch plot(fftshift(N)) 
y = (1000/2)*ifft(N);
subplot(2,1,2);
plot(y);

Rückwärts gerechnet. Wir wissen die halbe Amplitude entfällt auf jeweils 
eine Seite. Zudem kommt bei der diskreten Fouriertransformation die 
Anzahl der Samples als Umrechnungsfaktor hinzu. Da bei 1000 Samples nun 
mal die hälfte dem symmetrischen Anteil gehört, folgt: 1000/2 = 500

Und siehe da alles stimmt.

Mal anders herum gedacht, besitzt eine Seite deines Frequenzvektors auch 
nur N=500 Werte, somit hat er auch nur das halbe Gewicht.

Gruß Alexander

[Alexander L. (lippi2000)]

Nochmal, dein Spektrum geht immer von -fs/2 bis +fs/2.

Bisher wurde immer nur der Bereich 0...+fs/2 betrachtet. Dies entspricht 
deinen Frequenzvektor f(0:N/2-1). Woher nun den negativen Bereich 
nehmen???

Wir wissen, dass bei abgetasteten Signalen das Spektrum mit fs 
periodisch ist. Das heißt der Spektralanteil -fs/2...0 wird jetzt im 
Vektor f(N/2:N) sichtbar. Dieser ist demnach gegenüber dem f(0:N/2-1) 
gespiegelt, weswegen auch der Befehl "fliplr" verwendet wurde.

Was passiert, wenn du den Teil f(N/2:N) weglässt. Du hast ein 
einseitiges Spektrum. Diese rücktransformiert in den Zeitbereich ergibt 
ein analytisches Signal (komplexes Signal). Dies ist zwar in einigen 
Anwendungen (Quadratur-Signale) sinnvoll, aber entspricht nicht deinem 
gewünschten Signal.

[H. D. (hierunddort)]

Hi Alexander,

danke für deine Antworten und deine Geduld.

Der Matlab-Code funktioniert,
-allerdings sehe ich hier Widersprüche zu deinen Aussagen:
> dann muss jede Spektralseite nur die halbe Amplitude besitzen.
-Wieso nimmst du nur die Hälfte des von deinem Signal bekanntem 
Spektrums? Die Werte f(501:1000) lässt du völlig außer acht.
-Ich verstehe nicht weshalb der Umrechnungsfaktor halbiert wird.

Hier der Matlabcode für dein Beispiel:

N = zeros(1,1000);     
%Kosinusignal mit 10 Perioden in 1000 Samples 
%Achtung N(1) ist der DC-Wert
N(11) = 1;
%Es scheint auch ohne symmetrischen Teil zu funktionieren
subplot(2,2,1);
plot(N);          %Hier geht auch plot(fftshift(N)) 
y = (1000)*ifft(N);
subplot(2,2,3);
%Ich muss dazu lediglich den Realteil des Signals y auswerten
plot(real(y));

%mit symmetrischem Teil
N = [N(1) 0.5*N(2:500) 0.5*conj(fliplr(N(1:500)))];
subplot(2,2,2);
plot(N);          %Hier geht auch plot(fftshift(N)) 
y = (1000)*ifft(N);
subplot(2,2,4);
plot((y));

[H. D. (hierunddort)]

Meine Baustelle ist halt, dass ein Signal im Frequenzbereich vorgegeben 
ist. Sagen wir von f(0) bis f(500).
Ob nun f(500)=Fs oder f(500)=Fs/2 oder ganz was anderes ist nicht 
vorgegeben.

Um die meiste Information auszunutzen setze ich f(500)=Fs/2.
Demnach muss ich das Spektrum ins Zweiseitenband erweitern.
Wegen der Energieerhaltung, da dann jede Frequenz (außer f(0)) zweimal 
vorhanden ist, halbiere ich entsprechend
f_zsb=[f(1) 0.5*f(2:501) 0.5*conj(fliplr(f(2:500)))];

N = ones(1,500);     
%ein anderes Beispiel Signal
%N = N.*fliplr(1:500)/500;

%Es scheint auch ohne symmetrischen Teil zu funktionieren
subplot(2,2,1);
plot(N);axis([-100 1000 -0.1 1.1]); grid on;%
%ifft der Länge des Signals
y = ifft(N,500);
%Hier ist kein Umrechnungsfaktor nötig!!! 
subplot(2,2,3); 
%Ich muss dazu lediglich den Realteil des Signals y auswerten
stem(real(y)); axis([-100 1000 -0.1 1.1]); grid on;

%mit symmetrischem Teil
N = [N(1) 0.5*N(2:500) 0.5*conj(fliplr(N(2:500)))];
subplot(2,2,2); 
plot(N); axis([-100 1100 -0.1 1.1]); grid on;        %
%Hier ist kein Umrechnungsfaktor nötig!!! 
y = ifft(N,998);
subplot(2,2,4); 
stem(real(y)); axis([-100 1000 -0.1 1.1]); grid on;

Ich hoffe wir sind uns zumindest Einig, dass Spektrum und Signal von 
subplot 1 und subplot 3 übereinstimmen.

Weshalb fallen denn hier eigentlich die Umrechnungsfaktoren weg?

Demnach liegt irgendwo ein Fehler bei der Berechnungsmethode von Subplot 
2 und 4.

Dies wird noch deutlicher, wenn das zweite Beispiel Spektrum betrachtet 
wird, denn da ist es komplizierter.

Grüße Hierunddort

[Alexander L. (lippi2000)]

Überlege doch mal was bei der diskreten Fourier-Rücktransformation 
gemacht wird...

Die diskrete Fouriertransformation setzt ein abgetastetes Signal 
vorraus. Soweit klar.

Ein diskretes Signal besitzt allerdings ein periodisches Spektrum!!!

Wenn du jetzt die Rücktransformation durchführst, dann wird 
vorrausgesetzt, dass dein Frequenzfenster ebend periodisch ist.

Nun setzte deine beiden Spektren mal periodisch fort, indem du das 
Fenster von N jeweils einmal links und einmal rechts wieder anfügst.

Komisch, irgendwie sieht dein periodisches Spektrum beider Signale 
unterschiedlich aus.

Im ersten Fall: ähnelt einem periodischen "Sägezahn"
Im zweiten Fall: ähnelt einem periodischen "Dreieck"

Logisch ergeben sich auch zwei andere Signale.

Jetzt erweitere mal dein Plots um den Imaginärteil von y, dann wirst du 
sehen, dass sich im 1.Fall ein komplexes Signal ergibt und im 2.Fall ein 
reelles (ausgenommen von Rundungsfehlern).

Die Signalamplitude des 2.Signals ist nur halb so groß, da hier N=998 
Frequenzwerte rücktransformiert wurden anstatt 500.

[H. D. (hierunddort)]

Irgendwie hab ich das Gefühl wir Reden an einander vorbei...

ausgehend von:fliplr(1:500)/500
die Realteile von

ifft([fliplr(1:500)/500 zeros(1,500)])
und 
ifft([1 0.5*fliplr(2:500)/500 0.5*(1:500)/500])

stimmen überein, ganz gleich ob die Beträge der Spektren gleich 
aussehen, oder nicht.

Nach der zweiten Variante sind die Imaginärteile gering bis nicht 
vorhanden. Demnach wohl auch vorzuziehen.

Zurück zur Ursprünglichen Frage...
Ich habe einen einzelnen Peak X im Frequenzbereich.
um hiervon das Signal im Zeitbereich zu berechnen soll ich das Signal 
nach

x=ifft([1 0.5*X(1:499) 0.5*fliplr(conj(X(1:500)))])

berechnen und letztendlich mit der gesamten FFT-Länge multiplizieren.
Das klappt.


Noch eine (hoffentlich) abschließende Frage:
Wieso führt bei

X=ones(1,500);
X=[1 0.5*X(1:499) 0.5*fliplr(conj(X))]; % Das Signal entsprechend erweitert
x=ifft(X,1000);          % Reelles zeit Signal
x=x*1000;                %Umrechnung

die Multiplikation mit der gesamten FFT-Länge zu einer Amplitude die 
größer ist als ursprünglich erwartet?

Grüße Hierunddort

[Alexander L. (lippi2000)]

Also ich konnte bei beiden Varianten kein unterschied sehen.

Lediglich bein Variante 2 wird durch den Aufruf ifft(x,1000) bereits das 
Signal mit der Länge 1000 multipliziert.

Ich erhalte jeweils einen Dirac an der Stelle Null mit der Amplitude 
500.

[H. D. (hierunddort)]

Alexander Liebhold schrieb:
> Also ich konnte bei beiden Varianten kein unterschied sehen.
Sorry, jetzt bin ich mir nicht im Klaren darüber welche du meinst.
Falls

ifft([fliplr(1:500)/500 zeros(1,500)])
und 
ifft([1 0.5*fliplr(2:500)/500 0.5*(1:500)/500])

dann bitte nochmal nachrechnen

> Lediglich bein Variante 2 wird durch den Aufruf ifft(x,1000) bereits das
> Signal mit der Länge 1000 multipliziert.
Der Aufruf ifft(x,1000) legt doch lediglich fest, dass die ifft eine 
Länge von 1000 hat, zu kurze Signale werden mit 0 aufgefüllt, zu lanfe 
abgeschnitten.

> Ich erhalte jeweils einen Dirac an der Stelle Null mit der Amplitude
> 500.

Genau, Wieso 500?
Bei einem Sinus der Amplitude 1 im Zeitbereich, erhalte ich im 
Spektralbereich ein Zweiseitenband mit der Amplitude 0.5 (bzw. daraus 
ein Einseitenband der Amplitude 1)
Wieso ist bei einem Dirac-Impuls der Amplitude 1 im Zeitbereich, dann 
das Spektrum nicht 1 ?