WICHTIG: Der Betrieb von goMatlab.de wird privat finanziert fortgesetzt. - Mehr Infos...

Mein MATLAB Forum - goMatlab.de

Mein MATLAB Forum

 
Gast > Registrieren       Autologin?   

Partner:




Forum
      Option
[Erweitert]
  • Diese Seite per Mail weiterempfehlen
     


Gehe zu:  
Neues Thema eröffnen Neue Antwort erstellen

Wie FFT aus diskreten Datenpunkten (sinusähnlich) bestimmen

 

L. J.
Forum-Newbie

Forum-Newbie


Beiträge: 6
Anmeldedatum: 28.06.13
Wohnort: ---
Version: ---
     Beitrag Verfasst am: 30.07.2013, 01:52     Titel: Wie FFT aus diskreten Datenpunkten (sinusähnlich) bestimmen
  Antworten mit Zitat      
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
               


Lieder hat das mit der Interpolation nichts gebracht, wie man auf dem Bilder erkennen kann. Wenn ich nur
Code:
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?

Vielen Dank im Voraus
Gruß
L J


EDIT: Ohne es jetzt ganz verstanden zu haben, habe ich mal dieses Beispiel ausprobiert:
http://www.gomatlab.de/frequenz-aus.....t27707,highlight,fft.html
Das Ergebnis sieht dann so 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.
 
-> zweites Bild
Könnte mir jemand bitte erklären, was es mit dem NFFT auf sich hat?

fft2.png
 Beschreibung:

Download
 Dateiname:  fft2.png
 Dateigröße:  3.14 KB
 Heruntergeladen:  850 mal
fft.png
 Beschreibung:

Download
 Dateiname:  fft.png
 Dateigröße:  3.2 KB
 Heruntergeladen:  882 mal
Private Nachricht senden Benutzer-Profile anzeigen


Mmmartina
Forum-Meister

Forum-Meister


Beiträge: 745
Anmeldedatum: 30.10.12
Wohnort: hier
Version: R2020a
     Beitrag Verfasst am: 30.07.2013, 08:22     Titel:
  Antworten mit Zitat      
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))
Private Nachricht senden Benutzer-Profile anzeigen
 
L. J.
Themenstarter

Forum-Newbie

Forum-Newbie


Beiträge: 6
Anmeldedatum: 28.06.13
Wohnort: ---
Version: ---
     Beitrag Verfasst am: 30.07.2013, 11:31     Titel:
  Antworten mit Zitat      
Hallo Martina,

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.

plot(xfft)


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?

Gruß
L J
Private Nachricht senden Benutzer-Profile anzeigen
 
DSP
Forum-Meister

Forum-Meister



Beiträge: 2.117
Anmeldedatum: 28.02.11
Wohnort: ---
Version: R2014b
     Beitrag Verfasst am: 30.07.2013, 12:04     Titel:
  Antworten mit Zitat      
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?
if mod(nfft,2) == 0;
    k = (nfft/2) + 1;
else
    nfft = nfft + 1;
    k = (nfft/2) + 1;
end

fn = fa/2; % Nyquistfreq.
df = fa/nfft; % Frequenzauflösung des Spektrums
% Frequenzvektor: Darstellung bis Nyquistfreq.
fv = 0: df : fn;

sig = signal(:);

% Signal transformieren
Fy = fft(sig,nfft);
 
%   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)
 


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.
Private Nachricht senden Benutzer-Profile anzeigen
 
L. J.
Themenstarter

Forum-Newbie

Forum-Newbie


Beiträge: 6
Anmeldedatum: 28.06.13
Wohnort: ---
Version: ---
     Beitrag Verfasst am: 30.07.2013, 15:20     Titel:
  Antworten mit Zitat      
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([0 0.5]))


und die Funktion:
Code:
  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?
        if mod(nfft,2) == 0;
            k = (nfft/2) + 1;
        else
            nfft = nfft + 1;
            k = (nfft/2) + 1;
        end
   
    fn = fa/2; % Nyquistfreq.
    df = fa/nfft; % Frequenzauflösung des Spektrums
    % Frequenzvektor: Darstellung bis Nyquistfreq.
    fv = 0: df : fn;

    sig = signal(:);

    % Signal transformieren
    Fy = fft(sig,nfft);

    %   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.

fft3.png
 Beschreibung:

Download
 Dateiname:  fft3.png
 Dateigröße:  3.61 KB
 Heruntergeladen:  784 mal
Private Nachricht senden Benutzer-Profile anzeigen
 
DSP
Forum-Meister

Forum-Meister



Beiträge: 2.117
Anmeldedatum: 28.02.11
Wohnort: ---
Version: R2014b
     Beitrag Verfasst am: 30.07.2013, 15:36     Titel:
  Antworten mit Zitat      
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.
Private Nachricht senden Benutzer-Profile anzeigen
 
L. J.
Themenstarter

Forum-Newbie

Forum-Newbie


Beiträge: 6
Anmeldedatum: 28.06.13
Wohnort: ---
Version: ---
     Beitrag Verfasst am: 30.07.2013, 15:43     Titel:
  Antworten mit Zitat      
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?

Vielen Dank
L J

final.png
 Beschreibung:

Download
 Dateiname:  final.png
 Dateigröße:  81.03 KB
 Heruntergeladen:  758 mal
Private Nachricht senden Benutzer-Profile anzeigen
 
DSP
Forum-Meister

Forum-Meister



Beiträge: 2.117
Anmeldedatum: 28.02.11
Wohnort: ---
Version: R2014b
     Beitrag Verfasst am: 30.07.2013, 16:22     Titel:
  Antworten mit Zitat      
Zitat:

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.
Private Nachricht senden Benutzer-Profile anzeigen
 
L. J.
Themenstarter

Forum-Newbie

Forum-Newbie


Beiträge: 6
Anmeldedatum: 28.06.13
Wohnort: ---
Version: ---
     Beitrag Verfasst am: 30.07.2013, 23:14     Titel:
  Antworten mit Zitat      
danke für die Erklärung
Private Nachricht senden Benutzer-Profile anzeigen
 
Neues Thema eröffnen Neue Antwort erstellen



Einstellungen und Berechtigungen
Beiträge der letzten Zeit anzeigen:

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
.





 Impressum  | Nutzungsbedingungen  | Datenschutz | FAQ | goMatlab RSS Button RSS

Hosted by:


Copyright © 2007 - 2025 goMatlab.de | Dies ist keine offizielle Website der Firma The Mathworks

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.