Inverse Fouriertransformation in Matlab

Wissenschaft Mathematik
#1

Hallo,
stehe (mangels 100% Verständnis) auf dem Schlauch, beim Versuch die Matlab-Befehle für eine inverse Fouriertransformation zu schreiben.
Zur Übung wollte ich das Beispiel für eine DFT im Buch „Matlab-Simulink…“ von Angermann, etc. invertieren, also quasi bei gegebenem Spektrum das Signal anzeigen lassen. Allerdings macht mir die Rücktransformation der x-Achse zu schaffen.

t=0.01:0.01:0.5; %Zeitvektor
x=5+8*sin(2*pi*8*t)+4*cos(2*pi*33*t); %Signal
Ts=diff(t(1:2)); %Abtastzeit
N=length(x); %Länge des Datenvektors
f=[0:floor((N-1)/2)]/(N*Ts);
X=fft(x);
X=X/N;
X=[X(1) 2*X(2:floor((N-1)/2)+1)]; %Begrenzen auf

#2

Hallo,

zuerst fällt mir auf, dass du das Spektrum für die Abbildung veränderst. Erstens schneidest du die Hälfte ab da dich anscheinend nur die Amplitude interessiert. Diese Darstellung ist bei 1-dim Signalen üblich, aber es geht dadurch die Phase der negativen Freq. verloren, die in dem plot nicht erkennbar ist (denn die ist antisymmetrisch). Wenn du daher dieses halbierte Spektrum zurück transformierst, kommt daher was anderes raus. Da du weiter das Spektrum durch N teilst (warum eigtl?), würde ausserdem die Skalierung der rücktransformierten Originalfunktion verändert werden.
Die Freq.achse würd ich im plot einfach per Hand einstellen, und die Zeitachse einfach vom Original übernehmen. Die Zeit ist irrelevant für die Transformation. Wenn der betrachtete Zeitintervall die Länge T hat, dann ist die höchste Frequenz die du bekommen kannst 1/(2*T). T ist in meinem code range(t). Hoffe das hilft dir für deine Zwecke weiter. beste Grüsse, Florian

t=0.01:0.01:0.5; %Zeitvektor
x=5+8*sin(2*pi*8*t)+4*cos(2*pi*33*t); %Signal
N=length(x); %Länge des Datenvektors
f=linspace(-N/2, N/2, N)*1/(2*range(t)); %Zeitfrequenz
X=fft(x); %Signalspektrum

figure
subplot(131)
plot (t,x,’.-’)
subplot(132)
stem (f,fftshift(abs(X))) %wenn du hier nur das Spektrum für pos. Freq. zeigen willst, würde ich X in einer extra Variablen schneiden…
xlim([-N/2 N/2]);
subplot(133)
plot(t,ifft(X), ‚.-‘) %…denn sonst klappt die Rücktransformation nicht.

#3

Hi Florian,

vielen Dank für Deine schnelle Antwort.
Warum das Beispiel im Buch das Spektrum halbiert, weiss ich ehrlich gesagt auch nicht. Aber ich werde auf alle Fälle Deinen Code heute abend ausprobieren.

Mein eigenes Problem ist eigentlich, dass ich ein Spektrum S(w)=A*w^B und zugehörige Kreisfrequenzen w (von 0 bis 3,5 rad/sec) gegeben habe und jetzt gerne das zugehörige Signal ausrechnen und anzeigen würde. Danach will einen Filter auf das Spektrum legen und wieder das Signal anschauen - zwecks Feststellung von Unterschieden. Aber dank Deiner Hilfe komme ich jetzt hoffentlich weiter !

Besten Gruß
Peter

#4

Hi,

kann dir da leider „nicht mehr“ weiterhelfen, weil ich seit 3 oder 4 jahren nich mehr mit matlab arbeite.
Sorry.

gruß

#5

Ahoi,

aus einem gegebenen Spektrum bekommt man ganz leicht das Signal, indem man folgende Funktion nutzt:

f(t)=Summe(von k=0 bis N)(a(k)*cos(2*pi*f(k)*t)-b(k)*sin(2*pi*f(k)*t))

Dabei ist a(k)=real(X(k)) und b(k)=imag(X(k))

In Matlab entspräche das also:

u(1,1:length(t))=0; %Definition des Vektors u
for k=1:length(X)
u=u+real(X(k))*cos(2*pi*f(k)*t) -imag(X(k))*sin(2*pi*f(k)*t);
%u=u+real(X(k)*exp(2*pi*f(k)*1i*t)); %alternativ
end

figure
plot(t,u)

Wenn du das Ganze mit der ifft-Funktion machen möchtest, brauchst du alle Werte von X (also vor der Begrenzung und so. D.h. du musst zurückrechnen, was dank der Symmetrie möglich ist, oder in diesem Fall nimmst du einfach das X aus der FFT - also x=ifft(fft(x)) … naja meine Methode ist besser :smile:

Ich hoffe, ich konnte dir etwas helfen. Falls noch Fragen bestehen, frag ruhig - ich müsste mich aber ggf. nochmal etwas einlesen in die Materie. Sehr gut beschrieben sind die Fourier-Transformationen übrigens bei Wikipedia oder im Bronstein.

Gruß, Ulfson

#6

Hi,
Tut mir leid, aber ich hab Matlab das letzte mal vor 2 Jahren genutzt … Ich kann dir da leider nicht helfen.

MfG

#7

sorry bin schon länger raus aus dem Thema

#8

Hi Ulfson,

vielen Dank für Deine Hilfe - habe ich gerade ausprobiert und sieht gar nicht so schlecht aus !!!

Merci,
Flitzpiepe1

#9

@absOfZero

Hi,
jetzt bin ichs nochmal, da ich zu doof bin…
Also - ich habe einen Spektrumsverlauf S(w)(also kontinuierlich, noch nicht diskret) gegeben (w=Kreisfrequenz). Spektrum ist einseitig, sprich nicht symmetrisch zur y-Achse oder so… Ferner habe ich meinen gewünschten Kreisfrequenzrahmen von 0 bis 3.5 rad/sec (Schrittweite 0,00625 rad/sec).
Jetzt würde ich gerne a) das Spektrum (mit der Frequenz als x-Achse) anzeigen und b) das Signal anzeigen (mit der Zeit als x-Achse).
Wie schaut da mein Code aus bzw. wie transformiere ich die Achsen usw. …!?

Sorry - steh auf dem Schlauch…

Gruß
Peter

#10

Die Amplitude ist gleich für neg und pos Frequenzen (symmetrisch), die Phase ist negiert für negative Frequenzen (antisymmetrisch). Matlab verwendet nicht die Polarnotation (Amplitude und Phase), sondern Realteil und Imaginärteil (z.B. 4+2i). Hier gilt das selbe: Realteil symmetrisch, Imaginärteil asymmetrisch (auch wenn Realteil ungleich Amplitude und Imaginärteil ungleich Phase). Um das Spektrum zu rekonstruieren, musst du also den vorhandenen Teil für pos Frequenzen kopieren, spiegeln, seinen Imaginärteil negieren und an das ursprünglich vorhandene Spektrum dranhängen. Das ist bloss Gefummel. Die Umechnung Zeitachse-Frequenzachse musst du per Hand machen. Hier ein kleines Beispiel.

n = 100;
t = linspace(0,1,n); %Zeitachse
y = rand(1,n); %ein Rauschsignal

y_ = fft(y); %dessen Spektrum
f = linspace(-1/(2*range(t)),1/(2*range(t)), n+1); %die Frequenzachse - Min und Max sind gegeben durch die höchstmögliche Abtastfrequenz
f = f(2:n+1); %die höchste absolute Frequenz wird wie die Null nur einmal angegeben, ich leg sie mal in den pos Bereich

y_pos = y_(1:n/2+1); %das einseitige positive Spektrum
y_tmp = y_pos(n/2:-1:2); %das positive Spektrum gespiegelt ohne höchste und niedigste Frequenz
y_neg = real(y_tmp) - imag(y_tmp)*i; %Negieren des Imaginärteils
y_rec = [y_pos y_neg]; %Zusammensetzen - in Matlab werden die negativen Frequenzen an die positiven drangehängt

y_-y_rec %Testausgabe: sollten alles Nullen sein, wenn die Spektren identisch sind

figure
subplot(221)
plot(t,y);
xlabel(‚t‘);
ylabel(‚y‘);
title(‚Ausgangssignal‘);

subplot(222)
plot(f(find(f>=0)),abs(y_pos), ‚b‘);
xlim([0 max(f)]);
xlabel(‚f‘);
ylabel(‚Amplitude‘);
title(‚Original positives Spektrum‘);

subplot(223)
plot(t,ifft(y_rec));
xlabel(‚t‘);
ylabel(‚y‘);
title(‚Rekonstruiertes Signal‘);

subplot(224)
plot(f,abs(y_rec([n/2+2:n 1:n/2+1])));
xlabel(‚f‘);
ylabel(‚Amplitude‘);
title(‚Rekonstruiertes Spektrum‘);

Kann man bestimmt auch eleganter machen. Hoffe das hilft trotzdem erstmal weiter,
beste Grüsse
Florian

#11

Ok. - super Beispiel, so langsam wirds bei mir.
Aber da ich trotz allem immer noch Probleme habe, dreisterweise hier mein Code:

w=0:0.00625:3.5; %gegebene Kreisfrequenzen des Spektrums
n=length(w); %Anzahl Kreisfrequenzen
S=8.1e-3.*9.81^2./w.^5.*exp(-0.74.*((9.81/15.4./w).^4)); % mein einseitiges Spektrum
Stmp=S(n:-1:2); %Spektrum spiegeln
Srec=[S Stmp]; %Spektrum zusammensetzen
figure
subplot(211)
t=(-(n-1):1:frowning:n-1)); %Zeitachse
plot(t,ifft(Srec)) % Darstellung meiner Zeitfunktion
subplot (212)
plot(w,S) %Darstellung meines Spektrums

Leider bekomme ich dabei keine Zeitfunktion dargestellt !!???

Hoffentlich nicht nervend und mit bestem Dank
Peter

#12

Hallo,
Wenn ich nicht irre, ergibt sich der range der Zeitachse durch die niedrigste Frequenz. Die niedrigste (absolute) Frequenz schwingt genau einmal über den betrachteten Zeitintervall. Das ist hier wMin=0.00625 (die Frequenz Null gibt nur die baseline des Signals an, ihre Periodenlänge geht gegen unendlich). Damit kann man die obere Grenze des Intervalls berechnen: tMax=2*pi/wMin. Mit der unteren Grenze tMin=0 kann man dann mit linspace einen Zeitarray aufspannen.
Allerdings gibts hier ein paar Probleme: 1.) w=0 ist gegeben, das Spektrum ist für w=0 aber nicht gegeben da durch w geteilt wird. Man kann es als S(0)=0 „ergänzen“, das legt der plot jdf nahe. Die Form der Signalfunktion ändert sich durch die Frequenz 0 ohnehin nicht, sondern wird nur vertikal verschoben.
2.) Ich weiss nicht woher diese Aufgabe stammt, vielleicht beschreibt diese Funktion aber nur das Amlitudenspektrum, z.B. weil es sich um einen Filter (plot sieht so aus) handelt. In diesem Fall würde die Rückrechnung keinen Sinn machen, ausserdem ist es nicht möglich aus dem Amplitudenspektrum die Impulsfunktion zurückzugewinnen da die Information über die Phase fehlt. Der Verdacht liegt nahe, da die Spektrumsfunktion keinen Imaginärteil besitzt, das ist eher ungewöhnlich, ausserdem sieht die Impulsfunktion seltsam aus.

beste Grüsse
Florian

w = 0.00625:0.00625:3.5; %gegebene Kreisfrequenzen des Spektrums
S = [0 8.1e-3 * 9.81^2./w.^5.*exp(-0.74*((9.81/15.4./w).^4))]; % mein einseitiges Spektrum
n = length(S); %Anzahl Kreisfrequenzen
w = [0 w];
Stmp = S(n-1:-1:2); %Spektrum spiegeln
Srec = [S Stmp]; %Spektrum zusammensetzen

figure
subplot(211)
t=linspace(0,2*pi./min(w(find(w~=0))),2*n-2); %Zeitachse
plot(t,ifft(Srec)) % Darstellung meiner Zeitfunktion
subplot (212)
plot(w,S) %Darstellung meines Spektrums

#13

Hi,

also erstmal merci für Deine große Mühe. Wiklich sehr nett !
Ist mir zwar ein bisschen peinlich, aber inzwischen (und auch weil Du Deine Zweifel geäußerst hast) ist mir aufgegangen, dass es sich bei meinem Spektrum um das Leistungsdichtespektrum handelt !
Vielleicht muss ich dann erst das „normale“ Spektrum rückrechnen und mache dann mit meinem normalem Code weiter !?

Gruß
Peter

#14

Hallo,
laut Wikipedia ist das Leistungsdichtespektrum, hier auch Autoleistungsspektrum genannt, nicht gleich dem Frequenzspektrum des Signals, sondern ist dessen Power, und das ist wiederum das Quadrat des Betrages des Signal (d.h. das Quadrat der Amplitude): http://de.wikipedia.org/wiki/Spektrale_Leistungsdich…
Das bedeutet, dass die Phaseninformation nicht mit drin ist. Daher gibt es auch keinen imaginären Anteil im Spektrum. Da keine Information über die Phase vorliegt, ist es nicht möglich das Ausgangssignal zu rekonstruieren. Das ist wie wenn man versucht (a+b)^2=c zu lösen, wenn nur c gegeben ist; es gibt unendlich viele Lösungen.

beste Grüsse
Florian

1 Like
#15

Hi Florian

o.k. - gebe mich geschlagen ! Vielen Dank nochmal für Deine Hilfe. Hat mir echt sehr weiter geholfen, auch wenn ich ein bisschen langsam war !!!

Danke und Gruß
Peter