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:
1
t = 0:.001:.239;
2
punkte = length(t);
3
x = sin(2*pi*50*t);
4
subplot(3,1,1)
5
plot(t,x);
6
grid on;
7
8
points = 1000; % muss gerade und in N sein
9
10
%FFT
11
Y = fft(x,points)/punkte;
12
% Pyy=Y.*conj(Y)/points; %Leistungsdichtespektrum
13
f=1000*(1:(points/2))/points;
14
subplot(3,1,2)
15
plot(f,abs(Y(1:(points/2))));
16
grid on;
17
18
%Rücktransformation
19
f_p = [f fliplr(f)]; % Frequenzvektor
20
Y1 = Y.*f_p*2*pi*i;
21
T=ifft(Y1,points);
22
length(Y)
23
subplot(3,1,3)
24
plot((1:punkte)*0.001,T(1:punkte));
25
grid on;
Falls wer Vorschläge hat, bitte posten!
mfg,
Christoph
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
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
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_diff/fourier_diff.html
Das müsste man jetzt numerisch nachvollziehen.
mfg
Danke für die Mühe!
Aber irgendwie geht es noch immer nicht!
Code für Plot 1:
1
t = 0:.001:.239;
2
punkte = length(t);
3
x = sin(2*pi*40*t);
4
subplot(3,1,1)
5
plot(t,x);
6
grid on;
7
8
points = 10000; % muss gerade und in N sein
9
10
%FFT
11
Y = fft(x,points);
12
% Pyy=Y.*conj(Y)/points; %Leistungsdichtespektrum
13
f=1000*(1:(points/2))/points;
14
subplot(3,1,2)
15
plot(f,abs(Y(1:(points/2))));
16
grid on;
17
18
%Rücktransformation
19
f_p = [f fliplr(f)]; % Frequenzvektor
20
Y1 = Y;
21
T=ifft(Y1,points);
22
length(Y)
23
subplot(3,1,3)
24
plot((1:punkte)*0.001,real(T(1:punkte)));
25
grid on;
Dieses klappt einwandfrei!
Hier dann mit dem Ableiten:
Nur mehr der Code welcher anderes ist:
1
%Rücktransformation
2
f_p = [f fliplr(f)]; % Frequenzvektor
3
Y1 = Y.*f_p*2*pi*i;
4
T=ifft(Y1,points);
5
length(Y)
6
subplot(3,1,3)
7
plot((1:punkte)*0.001,real(T(1:punkte)));
8
grid on;
und dann noch für den Plot 3:
1
%Rücktransformation
2
f_p = [f fliplr(f)]; % Frequenzvektor
3
Y1 = Y.*f_p*2*pi*i;
4
T=ifft(Y1,points);
5
length(Y)
6
subplot(3,1,3)
7
plot((1:punkte)*0.001,imag(T(1:punkte)));
8
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
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...
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_diff_numeric/diskrete_fourier_diff_versuch2.html
mfg
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
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
Thread beobachten
Seitenaufteilung abschalten
