Skalierungsproblem? FFT eines Zeitsignals

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

Partner:

Forum

Option

[Erweitert]

• Diese Seite per Mail weiterempfehlen

Gehe zu:

Mike77

Forum-Anfänger


	Beiträge: 16

	Anmeldedatum: 03.12.08

	Wohnort: Berlin

	Version: 7.7.0 R2008b

Verfasst am: 25.08.2011, 19:33 Titel: Skalierungsproblem? FFT eines Zeitsignals

Hallo,

ich habe von einem Sinuszeitsignal eine FFT machen wollen.
Mir sind da gerade beim Vergleich MATLAB mit LabVIEW einige Verschiedenheiten aufgefallen, deshalb hab ich konkrete Frage zu den jeweiligen Plots die dieser Code erzeugt:

Code:

% Zeitbereich
% ----------------------------------

fa = 4096; % Abtastfrequenz
fn = fa/2; % Nyquistfrequenz
N = 2048; % gewünschte FFT-Länge (N=2^x, sonst wird der DFT-Algorithmus verwendet!)
df = fa/N; % Frequenzauflösung
% Erzeugung eines Datensatzes mit N Abtastwerten
% ----------------------------------------------
t = 0 : 1/fa : (N-1)/fa; % x-Vektor
% Frequenzvorgabe in Hz als ganzzahlig Vielfaches der Frequenzauflösung der DFT/FFT:
f1 = df*100; % bei fa = 8000 Hz und N = 1024 beträgt df = 7,8125 Hz und
% f1 damit 781,25 Hz
f1 = 784;
f1 = df;
phase = pi/2;

a1=2;
y = a1*sin(2*pi*100*t); % y-Vektor
% Graphische Darstellung
% ----------------------
% max. Amplitude zur Skalierung der graphischen Darstellung feststellen:
max_y = max(abs(y))*1.1;
fig = figure(1);
plot(y)
axis([-100 N+100 -max_y max_y])
title('Datensignal: Sinus mit 100 Hz, Amplitude = 2; N = 2048; f_a_b_t_a_s_t = 4096 Hz')
ylabel('Amplitude')
xlabel('Stützstelle N')
grid

% Frequenzbereich
% ----------------------------------

% Berechnung der FFT
% ------------------
H = fft(y, N);
% Berechnung des Amplitudengangs aus dem komplexen Frequenzvektor H:
amplH = abs(H);
% Amplitudenskalierung (Normierung auf N) und verschieben der Elemente des
% Amplitudenvektors, so dass die Darstellung des Amplitudengangs von -fn...0...fn
% erfolgen kann:
amplitudengang = fftshift(amplH/N);
% Graphische Darstellung
% ----------------------
% Frequenzvektoren (werden bei der graphischen Darstellung benötigt):
x_fn = 0 : df : fn-df;
x_fa = 0 : df : fa-df;
% max. Amplitude zur Skalierung der graphischen Darstellung feststellen:
%a = max([a1, a2, a3, a4, a5]); % wird später benötigt
a = a1;
fig = figure(fig+1);
stem(x_fa-fn, amplitudengang, 'b.-')
axis([-fn fn 0 a/2*1.1])
title('zweiseitiges FFT-Amplitudenspektrum')
ylabel('Amplitude')
xlabel(['Auflösung: ',num2str(df),' Hz, Frequenz in Hz'])
grid

% Darstellung des interessierenden Frequenzbereichs des
% Amplitudengangs (0...fn) und
% daran angepasste Amplitudenskalierung (Normierung auf N/2):
amplitudengang = [amplH(1)/N amplH(2:N/2)/(N/2)]; % DC-Bin auf N normieren!
fig = figure(fig+1);
stem(x_fn, amplitudengang, 'b.-')
axis([0 fn 0 a*1.1])
title('FFT-Spektrum / Betrag Amplitudengang (einseitig)')
ylabel('Amplitude')
xlabel(['Auflösung: ',num2str(df),' Hz, Frequenz in Hz'])
grid

% Ausgabe in dB
% ------------------
fig = figure(fig+1);
plot(x_fn, 20*log10(amplitudengang))
%axis([0 fn -100 20*log10(a)+3])
axis([0 2250 -400 30])
title('FFT-Spektrum / Betrag Amplitudengang (dB)')
ylabel('Amplitude in dB')
xlabel(['Auflösung: ',num2str(df),' Hz, Frequenz in Hz'])
grid

%Phasenspektrum
%-----------------------------------
fig = figure(fig+1);
H0 = fftshift(H);
f0 = (-N/2:N/2-1)*(fa/N);
p = angle(H0);
%p = unwrap(angle(H0));
f = 0 : df : 4094;
plot(f0,p*180/pi) % Ausgabe in Grad

title('FFT-Phasen-Spektrum')
xlabel('Frequenz (Hz)')
ylabel('Phase [°]')
grid

%Pwelch
%-----------------------------------

fig = figure(fig+1);
[Pxx,F]=pwelch(y,[],[],[],fa,'onesided'); % Uses default window, overlap & NFFT.
plot(F,10*log10(Pxx))
title ('Autoleistungsdichtespektrum (PSD)')
xlabel('Frequenz (Hz)')
ylabel('Leistungsdichte [dB/Hz]')
grid

%Leistungsdichte der DFT / FFT
%---------------------------------
lds=(abs(H).^2)/N;
% oder so: power = H.*conj(H)/N; ist das gleiche!

Leistungsspektrum = H.*conj(H)/N;

fig = figure(fig+1);
H0 = fftshift(H); % Rearrange H values
f0 = (-N/2:N/2-1)*(fa/N); % 0-centered frequency range
power0 = H0.*conj(H0)/N; % 0-centered power
plot(f0,power0)*(N/fa);
xlabel('Frequency (Hz)')
ylabel('Power')
title('{\bf 0-Centered Periodogram}')
grid

Funktion ohne Link?

Plot1: alles okay (auch in LabVIEW so)

Plot2: 2 Peak-Amplituden der höhe 1, passt auch

Plot3: jetzt einseitiges Betragsspektrum Amplitude also mal 2 = 1*2 = 2, passt auch. Hab ich auch so in LabVIEW)

Plot 4: Stell ich das nun als Dezibel dar ... da ist meine Spitze bei 6,02
bei LabVIEW sieht das generell gleich aus nur der peak ist bei 3,01 ?? hmm..und im unteren Bereich geht er bis -350, hier in MATLAB eher bis -320 dB?

Plot5: hier hab ich zweiseitig die Phase mal geplottet. Das ist sone sache, LanVIEW zeigt auch so eine Art Rauschen der Phase an. Springt immer zwischen 180 und -180 ca.? Wie muss das bei dem normalen Sinus mit 100 Hz aussehen? Weil wenn ich das mittel iwie, zeigt er mit bei der Phase in LabVIEW nichts mehr an. Also alles Null? Wie muss es korrekt sein?

Plot6: PSD, da weiss ich nicht mit den Parametern ist das so korrekt? LabVIEW geht da bei mir auch bis -320 dB.. hier ja nur bis -60. Hab Hamming fenster bei beiden und db! Oder liegt es an der Mittelungsmethode??

Plot7: das ist mir das größe Rätsel. Ich will einfach nur das Leistungsdichtespektrum. Das heisst ich nehm doch den Betrag des FFT-Wertes und quadriere ihn!? Dann bekomm ich also einen Peak der die Höhe 2 hat bei 100 Hz, oder? Zeigt mir LabVIEW auch so an. MATLAB macht hier iwas mit 2000?? 2048 iwie...sieht aus als hätte das mit meiner FFT Länge zu tun der faktor..

Kann jemand Tips geben?
Danke.

Atmos_kk

Forum-Fortgeschrittener


	Beiträge: 93

	Anmeldedatum: 23.05.11

	Wohnort: ---

	Version: ---

Verfasst am: 31.08.2011, 08:09 Titel:

Moin,
ich bin mir nicht ganz sicher was du vorhast. Was hast du für ein Signal zum verarbeiten?
Was ich aber sicher sagen kann ist, dass du für die dB-Umrechnung auf jeden Fall einen Referenzdruck benötigst. Der fehlt in dem FFT-Beispiel.

Greetz

DSP

Forum-Meister


	Beiträge: 2.117

	Anmeldedatum: 28.02.11

	Wohnort: ---

	Version: R2014b

Verfasst am: 31.08.2011, 09:31 Titel:

Bei Plot 4 kann ich dir zumindest helfen:

Es kommt einfach darauf an, um du Leistungs- oder Spannungspegel darstellen willst.

Bei Spannung: 20*log10(U/U0)
Leistung: 10*log10(P/P0)

Eine Amplitude von 2 ergibt daher ca. 6 bzw 3 dB je nach verwendeter dB Umrechung.

Bei dem Phasenspektrum kann ich keinen Fehler in der Berechung feststellen.


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 |

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.