Verfasst am: 30.07.2013, 01:52
Titel: Wie FFT aus diskreten Datenpunkten (sinusähnlich) bestimmen
Hallo,
ich komme gerade mit meiner Fouriertransformation nicht weiter.
Ausgangspunkt ist ein Array aus einzelnen Datenpunkten (ein Helligkeitsverlauf zwischen Hell und Dunkelflächen).
Die Datenpunkte ergeben geplottet ein sinus-ähnliches Signal, daher hatte ich gehofft darauf einfach eine FFT anwenden zu können.
Die direkte Anwendung der FFT auf das Signal brachte nichts brauchbares; zunächst vermutete ich, dass dies an den diskreten Werten liegen könnte und versuchte es dann mit einer Interpolation:
Code:
matrixwerte=[double(im_gray(ze, round(0.925*((right-left)/2)):round(1.075*((right-left)/2))));1:(abs(round(0.925*((right-left)/2))-round(1.075*((right-left)/2))))+1] %ergibt den Ausschnitt auf dem oberen Bild, hier mit zusätzlichen X-Koordinaten
y01=matrixwerte(1,:)/2^16; %y-werte
x01=matrixwerte(2,:); %x-werte
Interpolationslinie=pchip(x01, y01, 1:(abs(round(0.925*((right-left)/2))-round(1.075*((right-left)/2))))+1) %1. Ansatz: Interpolation satt diskreter Werte
eingebe, zeichnet er mir immerhin ein Dreieck;)
Es müsste doch irgendwie möglich sein, diese Modulation im Ortsraum in den Frequenzraum zu bekommen.
Wo ist hier meine Denkfehler (Fehlen Grenzen, muss ich noch irgendetwas definieren,...), oder welche Funktion könnte mir weiter helfen?
T0 = (x01(end)-x01(1))/(length(x01)-1); % sampling time
Fs = 1/T0; % ampling frequency
L = length(x01); % length of the signal
NFFT = 2^nextpow2(L); % next higher power of 2 of the length L (real) FFT = fft(y01,NFFT)/L; % FFT of x
freq = Fs/2*linspace(0,1,NFFT/2); % calculation of the frequency vector
xfft = 2*abs(FFT(1:NFFT/2)); % Calculate single-sided amplitude spectrum.
Was genau möchtest du mit der FFT erreichen?
Falls du aus der FFT die Frequenz/Amplitude, etc. herausrechen willst, ist das Möglich, aber nicht so trivial.
Eine recht gute und schnelle Methode, um dies für Sinussignale zu machen ist aber der Prony - allgorithmus (einfach mal google befragen, gibts sogar schon fertig als Matlab Code).
Zu deinen Ergebnissen:
Ich sehe da im ersten FFT-Bild eine ordentlichen FFT eines annähernden Sinussignals.
Beim Zweiten FFT-Bild hast du meiner Meinung nach die Achsen vertauscht - daher auch dieses komische gekruschel im Bild - wie sehen denn die Daten aus?
_________________
LG
Martina
"Wenn wir bedenken, daß wir alle verrückt sind, ist das Leben erklärt." (Mark Twain))
danke für die schnelle Antwort.
Ich sehe gerade, dass die Bilder vertauscht sind, also das erste Bild (fft2.png) ist das Ergebnis aus
Code:
T0 = (x01(end)-x01(1))/(length(x01)-1); % sampling time
Fs = 1/T0; % ampling frequency
L = length(x01); % length of the signal
NFFT = 2^nextpow2(L); % next higher power of 2 of the length L (real) FFT = fft(y01,NFFT)/L; % FFT of x
freq = Fs/2*linspace(0,1,NFFT/2); % calculation of the frequency vector
xfft = 2*abs(FFT(1:NFFT/2)); % Calculate single-sided amplitude spectrum.
Mir geht es eigentlich nur um die qualitative Darstellung der FT - mit den Werten der FT soll vorerst nichts weiter gemacht werden.
Das Ergebnis des o.g. Codes sieht ja an sich plausibel aus - wenn es soweit korrekt ist, würde es mir für die Darstellung bereits ausreichen. Lässt sich die Auflösung eventuell noch erhöhen und ist die x-Achse nun schon in f=1/px?
Die Amplitude ist mit dem Code so nicht ganz korrekt beim Gleichsignalanteil.
Code:
function[mag, mag_dB, fv] = FFT_betragsspektrum( signal, nfft, fa, scale) ;
% Input: % Signal im Zeitbereich % nfft = Anzahl Messwerte für fft % wenn nfft > length(sig) -> fft(sig,nfft) führt Zeropadding durch % fa = Abtastfreq. % scale: 0 = keine Impulsantwort als Eingang, 1 = Impulsantwort % Output: % Magnitude des Spektrums linear und dB skaliert % Frequenzvektor fv in [Hz] von 0...fa/2
% un-,gerade Anzahl Messwerte? ifmod(nfft,2) == 0;
k = (nfft/2) + 1;
else
nfft = nfft + 1;
k = (nfft/2) + 1;
end
Das Interpolieren bringt dir bei der FFT eigentlich nicht viel, da keine neuen Informationen entstehen und zusätzlich Aliasing Artefakte auftauchen werden. Man kann die Frequenzauflösung verbessern, in dem man die Anzahl nfft erhöht und dadurch automatisch Nullen beim FFT Befehl angefügt werden. Siehe Code...
Ich würde das Signal aber eher als ein Rechtecksignal denn ein Sinus sehen.
Danke für den Code.
So ganz verstehe ich die Schritte noch nicht, aber ich habe dennoch versucht das Ganze in meine Programm einzubinden.
Hier das Ergebnis:
Code:
% Berechnungen für FFT an Kante
matrixwerte=[double(im_gray(ze, round(0.925*((right-left)/2)):round(1.075*((right-left)/2))));1:(abs(round(0.925*((right-left)/2))-round(1.075*((right-left)/2))))+1];
y01=matrixwerte(1,:)/2^16; %y-werte
x01=matrixwerte(2,:); %x-werte
T0 = (x01(end)-x01(1))/(length(x01)-1); % sampling time
Fs = 1/T0; % sampling frequency Ls = length(x01); % length of the signal
NFFT = 2^nextpow2(Ls); % next higher power of 2 of the length L (real)
FTwerte= FFT_betragsspektrum( y01, NFFT, Fs, 0)
FTwerte(1,:)
FTwerte(2,:) subplot(4,4,8) %FT der Kante
plot(FTwerte(2,:),FTwerte(1,:)) title(' FT an steiler Kante' , 'FontWeight','bold');
xlabel('Frequenz [1/px]') xlim([00.5]))
function[plotwerte] = FFT_betragsspektrum( signal, nfft, fa, scale) % Input: % Signal im Zeitbereich % nfft = Anzahl Messwerte für fft % wenn nfft > length(sig) -> fft(sig,nfft) führt Zeropadding durch % fa = Abtastfreq. % scale: 0 = keine Impulsantwort als Eingang, 1 = Impulsantwort % Output: % Magnitude des Spektrums linear und dB skaliert % Frequenzvektor fv in [Hz] von 0...fa/2
% un-,gerade Anzahl Messwerte? ifmod(nfft,2) == 0;
k = (nfft/2) + 1;
else
nfft = nfft + 1;
k = (nfft/2) + 1;
end
% Betrag - nur positives Freq.spektrum if scale == 0
Fy_pos0 = abs(Fy(1:k));
% Skalierung
mag = [Fy_pos0(1)/nfft ;Fy_pos0(2:k-1)/(nfft/2);Fy_pos0(k)/nfft];
else% nicht durch nfft teilen bei Impulsantwort
mag = abs(Fy(1:k));
end
% Umrechnung in dB
mag_dB = 20*log10(mag + eps); % eps = kleine Konstante zur Vermeidung von log(0
plotwerte(1,:)=transpose(mag);
plotwerte(2,:)=fv;
end
Das Ergebnis sieht an sich nicht schlecht aus, nur etwas grob - kommen die Amplituden jetzt hin?
Man sieht eine "Hauptperiode" von etwa T=30px, also hätte ich bei f=1/30=0,03 etwas erwartet und dort liegt auch der zweite Peak. Entsteht der Wert bei 0 aufgrund des endlichen Signals->~Periode->infty; f->0?
Was das Eingangssignal angeht, hast du Recht, es ist eigentlich mehr ein Rechteck (so sind auch die Übergänge im Stern), aber bei schlechten Objektiven wird es schon eher etwas Sinus-ähnliches und wenn der Bildprozessor dann noch nachschärft wirds noch lustiger.
Wenn du eine höhere Auflösung willst, musst du NFFT größer machen. Die Abtastfrequenz kannst du ja schlecht ändern (außer druch Interpolation...aber das hatten wir schon).
Ein Rechtecksignal lässt sich doch mittels FFT mit vielen einzelnen Sinusschwingungen nachbilden. Die Grundfreq. deines Rechtecksignals liegt wohl bei ca. 0.04. Die anderen Anteile, davor und dannach geben dem Gleichsignalanteil bzw. höhere Frequenzanteile wieder mit denen das Dach "nachgebildet" wird.
Außerdem tritt hier mit Sicherheit auch der Leakage (Leck) Effekt auf, wodurch das Spektrum neben der groben Auflösung zusätzlich verschmiert wird. Das ließe sich nur durch eine vorherige Fensterung des Zeitsignals minimieren.
ok, vielen Dank.
Dann wirds wohl soweit passen. In der logarithmischen Darstellung erkennt man auch deutlich mehr.
Was hat es denn nun mit der NFFT auf sich - wieso wurde in dem ersten Beispiel 2^nextpow2(Ls) verwendet?
Noch eine kleine OT-Frage, Mein Auswerteskript wird für meherere Bilder durchlaufen, dabei entsteht jedes Mal eine Figure. Kann ich diese Figures am Ende irgendwie "einsammeln" und einer GUI mit Tabs darstellen?
Was hat es denn nun mit der NFFT auf sich - wieso wurde in dem ersten Beispiel 2^nextpow2(Ls) verwendet?
Damit der schnelle FFT Algorithmus verwendet wird, muss NFFT eine 2er Potenz sein. Ansonster wird der langsamere DFT Alg. verwendet. Ist NFFT > length(signal) werden automatisch die fehlenden Messwerte zur nächsten 2er Potenz durch Nullen angefügt.
Du kannst Beiträge in dieses Forum schreiben. Du kannst auf Beiträge in diesem Forum antworten. Du kannst deine Beiträge in diesem Forum nicht bearbeiten. Du kannst deine Beiträge in diesem Forum nicht löschen. Du kannst an Umfragen in diesem Forum nicht mitmachen. Du kannst Dateien in diesem Forum posten Du kannst Dateien in diesem Forum herunterladen
MATLAB, Simulink, Stateflow, Handle Graphics, Real-Time Workshop, SimBiology, SimHydraulics, SimEvents, and xPC TargetBox are registered trademarks and The MathWorks, the L-shaped membrane logo, and Embedded MATLAB are trademarks of The MathWorks, Inc.