mikrocontroller.net

Forum: Digitale Signalverarbeitung / DSP Matlab: fft() und ifft()


Autor: Christoph Gradl (christoph_gradl)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
Hallo,
kurz zu meinem Problem:
Ich möchte Messdaten (Positionsdaten) mit Hilfe von fft() 
transformieren. Im Frequenzbereich möchte ich dieses dann mal j*omega 
rechnen, welches ja einer Zeitableitung gleichkommt. Dieses Spektrum 
möchte ich wieder in den Zeitbereich zurück transformieren mit ifft(). 
Somit hätte ich so, ohne weiteres Filtern usw., die Geschwindigkeit.

Jedoch funktioniert dieses in meinem Fall nicht. Die Amplitude der 
"Geschwindigkeit" stimmt nicht überein.

Hier mein Code:
t = 0:.001:.239;
punkte = length(t);
x = sin(2*pi*50*t);
subplot(3,1,1)
plot(t,x);
grid on;

points = 1000; % muss gerade und in N sein

%FFT
Y = fft(x,points)/punkte;
% Pyy=Y.*conj(Y)/points; %Leistungsdichtespektrum
f=1000*(1:(points/2))/points;
subplot(3,1,2)
plot(f,abs(Y(1:(points/2))));
grid on;

%Rücktransformation
f_p = [f fliplr(f)]; % Frequenzvektor
Y1 = Y.*f_p*2*pi*i;
T=ifft(Y1,points);
length(Y)
subplot(3,1,3)
plot((1:punkte)*0.001,T(1:punkte));
grid on;

Falls wer Vorschläge hat, bitte posten!

mfg,
Christoph

Autor: R. M. (exp)
Datum:
Angehängte Dateien:

Bewertung
0 lesenswert
nicht lesenswert
Hallo,

wenn ich das richtig verstehe möchtest Du Dir das Differenzieren sparen.

Normalerweise verwendet man dazu aber nicht die FourierTransformation

sondern die LaplaceTransformation.

Wie meine angehängte Simulation aber zeigt, könnte man auch mit der FT

zum Ziel kommen, wenn von der Rücktransformierten Real- und Imaginärteil

vertauscht werden.

Ob das aber allgemeingültig ist, müsste man noch untersuchen,

ich will es keinesfalls behaupten.


P.S.
matlab habe ich leider nicht, ich simuliere alles mit MuPad.

mfg

Autor: Christoph Gradl (christoph_gradl)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
Hallo,
ich habe einen Vektor mit Messwerten und aus diesen möchte ich zeitlich 
ableiten. Die Laplacetransformation kann ich nur nehmen wenn ich 
symbolische Ausdrücke habe. Ist das richtig bzw. habe ich da etwas 
falsch in Erinnerung? Also bleibt mit nur mehr die fft() mit welcher ich 
dies lösen könnte.

In Maple habe ich die selben Plots wie du gemacht, darum habe ich mir 
gedacht das dies auch auf numerischer Basis vielleicht in Matlab möglich 
ist und ich könnte auf die Geschwindigkeit kommen ohne die Messdaten zu 
filtern.

mfg

Autor: R. M. (exp)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
Ja, Du hast schon recht, die LaplaceTransformation ist meines Wissens
auch nur für symbolische Ausdrücke anwendbar.
Jedenfalls kenne ich keinen numerischen Algo dafür.

Darum versuche ich immer zuerst mit meinem CAS eine Lösung zu finden,
die ich dann numerisch mit C++ oder Java implementieren kann ( mangels 
fehlendem MATLAB ).

In diesem Fall schauts wirklich so aus als ob folgendes zum Ziel führen 
würde:

X = FFT(x)
-> X' = X * s (multiplikation im Bildbereich mit s=jW)
-> Y = INVFFT(X')

und dann entspicht der IMAGINÄRTEIL(Y) der Ableitung von x.

Ich habe es mal mit ein paar Standart Funktionen ausprobiert
und die Ableitung wurde immer richtig berechnet

http://home.mnet-online.de/reimay/MuPad/fourier_di...

Das müsste man jetzt numerisch nachvollziehen.

mfg

Autor: Christoph Gradl (christoph_gradl)
Datum:
Angehängte Dateien:

Bewertung
0 lesenswert
nicht lesenswert
Danke für die Mühe!
Aber irgendwie geht es noch immer nicht!

Code für Plot 1:
t = 0:.001:.239;
punkte = length(t);
x = sin(2*pi*40*t);
subplot(3,1,1)
plot(t,x);
grid on;

points = 10000; % muss gerade und in N sein

%FFT
Y = fft(x,points);
% Pyy=Y.*conj(Y)/points; %Leistungsdichtespektrum
f=1000*(1:(points/2))/points;
subplot(3,1,2)
plot(f,abs(Y(1:(points/2))));
grid on;

%Rücktransformation
f_p = [f fliplr(f)]; % Frequenzvektor
Y1 = Y;
T=ifft(Y1,points);
length(Y)
subplot(3,1,3)
plot((1:punkte)*0.001,real(T(1:punkte)));
grid on;

Dieses klappt einwandfrei!

Hier dann mit dem Ableiten:
Nur mehr der Code welcher anderes ist:
%Rücktransformation
f_p = [f fliplr(f)]; % Frequenzvektor
Y1 = Y.*f_p*2*pi*i;
T=ifft(Y1,points);
length(Y)
subplot(3,1,3)
plot((1:punkte)*0.001,real(T(1:punkte)));
grid on;


und dann noch für den Plot 3:
%Rücktransformation
f_p = [f fliplr(f)]; % Frequenzvektor
Y1 = Y.*f_p*2*pi*i;
T=ifft(Y1,points);
length(Y)
subplot(3,1,3)
plot((1:punkte)*0.001,imag(T(1:punkte)));
grid on;

Warum muss man eigentlich auch den Imaginärteil nehmen?

Beschreibung der Plots:
1. Teil: Eingangssignal welches transformiert werden soll
2. Teil: Spektrum
3. Teil: Zeitsignal nach der Rücktransfromation


mfg

Autor: Flutsch (Gast)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
moin,
ohne das jetzt mal ausprobiert zu haben...könnte es vlt. sein das dein 
problem der Leck-Effekt ist?
weil zumindest in deinen bildchen sieht's so aus als ob dein Messlänge 
kein vielfaches der Periodendauer des Sinus sind und genau das sieht 
nach Leck-Effekt aus...daher auch im Spektrum nicht der reine 
dirac-impuls sondern die si-Fkt....
kann aber auch sein das ich mich irre =)

ich würd's einfach mal testen indem du nicht alle deine Messdaten 
verarbeitest, sondern die letzte unvollständige Periode einfach mal 
weglässt...wenn's damit schon gelöst sein sollte kannste ja immernoch 
entscheiden ob du wirklich alle Messdaten verarbeiten willst/musst oder 
nicht... bei ersterem müsstest du dich dann noch mit Fensterfunktionen 
auseinandersetzen...

Autor: Flutsch (Gast)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
also wenn ich's richtig seh währe es bei 0.225s sprich das neunfache der 
Periodendauer

Autor: R. M. (exp)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
Hallo
der 2. Plot sieht doch schon ganz gut aus.

Bis auf die Ein- und Ausschwinger ist der Ausgang gegenüber dem Eingang 
um 90 Grad verschoben. ( Wenn ich das richtig sehe )

Ich habe auch ein paar numerische Simulationen gemacht, es ist doch der
Realteil der ausgewertet werden soll.


http://home.mnet-online.de/reimay/MuPad/fourier_di...

mfg

Autor: Michael Lenz (hochbett)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
Hallo Christoph,

Dein Problem besteht darin, daß Du die positiven und die negativen 
Frequenzen nicht getrennt behandelst.

Die Fourierkoeffizienten der positiven Frequenzen mußt Du um 90° drehen, 
die der negativen Frequenzen um -90°.

------------------------------------------------------------------------
%% Datendefinition
winkel = 90;             % Drehwinkel in Grad
fs = 1500;                % Abtastfrequenz
f0 = 50;                 % Signalfrequenz
N  = 128;
t=0:1/fs:(N-1)/fs;       % Zeitachse
f  = 0:fs/N:(N-1)*fs/N;  % Frequenzachse

x1 = 1+sin(2*pi*f0*t);       % Sinussignal


%% Transformation
ft = fft(x1);

if rem(N,2) == 0,        % falls N gerade ist
    %ft(1) = 0;          % Gleichanteil, ggf. auskommentieren
    ft(2:round(N/2)+1)        = ft(2:round(N/2)+1) * 
exp(j*winkel/180*pi);
    ft(N:-1:round(N/2)+2) = ft(N:-1:round(N/2)+2) * 
exp(-j*winkel/180*pi);
else                     % falls N ungerade ist
    %ft(1)=0;            % Gleichanteil, ggf. auskommentieren
    %ft(ceil(N/2))=0;    % Gleichanteil, ggf. auskommentieren
    ft(2:ceil(N/2)-1)       = ft(2:ceil(N/2)-1) * exp(j*winkel/180*pi);
    ft(N:-1:ceil(N/2)+1)    = ft(N:-1:ceil(N/2)+1) * 
exp(-j*winkel/180*pi);
end;

x2 = ifft(ft);           % Rücktransformation

plot(t,x1,'b-'); hold on;
plot(t,real(x2),'r-'); hold on;                         % 
Realteilbildung wegen Rechenungenauigkeiten
plot(t,ft(1)/N-imag(hilbert(x1)),'g--'); hold off;      % Lösung über 
Hilberttransformation zum Vergleich

Gruß,
  Michael

Autor: Christoph Gradl (christoph_gradl)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
Hallo,
da habe ich wohl das Vorzeichen der negativen Frequenzen unter den Tisch 
fallen lassen. :-(

Vielen Dank für den Hinweis, schaut so aus als ob es jetzt funktioniert.

mfg,
Christoph

Antwort schreiben

Die Angabe einer E-Mail-Adresse ist freiwillig. Wenn Sie automatisch per E-Mail über Antworten auf Ihren Beitrag informiert werden möchten, melden Sie sich bitte an.

Wichtige Regeln - erst lesen, dann posten!

  • Groß- und Kleinschreibung verwenden
  • Längeren Sourcecode nicht im Text einfügen, sondern als Dateianhang

Formatierung (mehr Informationen...)

  • [c]C-Code[/c]
  • [avrasm]AVR-Assembler-Code[/avrasm]
  • [code]Code in anderen Sprachen, ASCII-Zeichnungen[/code]
  • [math]Formel in LaTeX-Syntax[/math]
  • [[Titel]] - Link zu Artikel
  • Verweis auf anderen Beitrag einfügen: Rechtsklick auf Beitragstitel,
    "Adresse kopieren", und in den Text einfügen




Bild automatisch verkleinern, falls nötig
Bitte das JPG-Format nur für Fotos und Scans verwenden!
Zeichnungen und Screenshots im PNG- oder
GIF-Format hochladen. Siehe Bildformate.
Hinweis: der ursprüngliche Beitrag ist mehr als 6 Monate alt.
Bitte hier nur auf die ursprüngliche Frage antworten,
für neue Fragen einen neuen Beitrag erstellen.

Mit dem Abschicken bestätigst du, die Nutzungsbedingungen anzuerkennen.