mikrocontroller.net

Forum: Digitale Signalverarbeitung / DSP Matlab Problem mit FFT


Autor: Zonck (Gast)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
Hallo Leute.
Ich mache gerade erste Gehversuche in Matlab und habe mit einen Sinus 
generiert, den ich gerne mittels FFT in den Spektralbereich 
tranformieren möchte. Das der Sinus eine Frequenz von 50Hz hat, erwarte 
ich natürlich auch einen Dirac bei 50Hz im Spektrum...aber den bekomme 
ich nicht - sondern bei etwa 54Hz! Warum?
Danke für euer Bemühen!!!

Hier mein Code:

t = 0:.001:.239;
x = sin(2*pi*50*t);
subplot(2,1,1)
hndl=plot(t,x)
grid on;

%FFT vom Signal x
Y= fft(x,240);
Pyy=Y.*conj(Y)/240; %Leistungsdichtespektrum
f=1000/240*(1:120);
subplot(2,1,2)
hndl=plot(f,Pyy(1:120));
grid on;

Autor: Thomas (Gast)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
Erhöhe mal die Genauigkeit, dann landet er auch bei der 50. Also von den 
240 z.B. auf 2400 gehen.

Autor: Zonck (Gast)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
Hallo nochmal, und danke Thomas!

Ich hatte die Genauigkeit schon etwas hochgedreht aber wie dich gezeigt 
hat, muss man sie "massiv" hochdrehen. Nun habe ich das Spektrum so, wie 
es sein soll. Doch bei der erneuten Rücktransformation (es sollte ja 
wieder das Zeitsignal rauskommen) klappt wieder was nicht - glaube 
Matlab will mich ärgern! :-(

t = 0:.001:.239;
x = sin(2*pi*50*t);
subplot(3,1,1)
hndl=plot(t,x)
grid on;

%FFT vom 10000);
Pyy=Y.*conj(Y)/10000; %Leistungsdichtespektrum
f=1000/10000*(1:5000);
subplot(3,1,2)
hndl=plot(f,Pyy(1:5000));
grid on;

%Rücktransformation
T=ifft(Y);
subplot(3,1,3)
hndl=plot(T);
grid on

Problem: Ich denke, die Zeitachse ist falsch skaliert!?

Danke

Autor: Thomas (Gast)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
Tach.

Du willst wohl folgendes:

t = 0:.001:.239;
punkte = length(t);
x = sin(2*pi*50*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,Pyy(1:(points/2)));
grid on;

%Rücktransformation
T=ifft(Y,points);
length(Y)
subplot(3,1,3)
plot(1:punkte,T(1:punkte));
grid on;


Beachte den Unterschied zwischen "punkte" und "points", ungeschickte 
Namenswahl, ok, bin schuld. Aber da siehst du die Skalierung.

Autor: Dennis (Gast)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
Hallo,

ich beschäftige mich momentan auch mit der FFT in Matlab.
Dabei sind nun drei Fragen aufgetreten.
Vielleicht kann mir jemand helfen.

1.Ich habe das Programm ähnlich aufgebaut. Nun verstehe ich aber zwei 
Zeilen nicht.
f=1000*(1:(points/2))/points;
Was passiert hier genau?
plot(f,Pyy(1:(points/2)));
hier wird ja der Graph erzeugt. Doch was bewirkt (1:(points/2)) ?

2. Die FFT ist ein Algorithmus der DFT. Bei der die Abtastrate eine 
zweierpotenz sein muss. Aber ist es ja in diesem Fall nicht, weil mit 
10000 abgetastet wird. ???

3. Mein Spektrum passt soweit. Was ich nicht verstehe ist, warum am Peak 
eine Abweichnung von +-1hz vorliegt. z.B. Frequenz ist auf 50Hz 
eingestellt.
Nun fängt das Spektrum bei 49Hz an zu steigen und hat bei 50Hz sein 
Maximum und fällt dann bis 51Hz.
Egal was ich änder, die +-1Hz bleiben. Wieso?

Kann mir da jemand weiterhelfen?

Autor: die ??? (Gast)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
1. Var(x:y) == Var ist ein Vector*, schneide den Bereich mit den Indizes
               x bis y aus.

 *) Kann auch ein mehrdimensionales Array sein.


2. Ist MATLAB egal, da wird zero-padding durchgeführt (Nullen anfügen
   bis zur nächsten 2er-Potenz).  help fft


3. Leckeffekt. Ausserdem verbindet MATLAB mit dem Aufruf plot alle
   Datenpunke mit einer Linie, was zu vermeindlichen 'Flanken' führt.
   Versuch mal stem.

Autor: Dennis (Gast)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
Alles klar.

Danke!

Autor: Dennis (Gast)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
Hallo,
ich habe ein Verständnis Problem mit meiner Ausgabe. Ich erhalte zwei
Spektren, wenn ich eine Sinus Funktion mit der FFT transformiere.
Vielleicht kann mir jemand erklären, wieso ich ein weiteres Spektrum bei
900Hz erhalte. (Bzw. bei Abtast-f)
Hier mein Programm:

clc

%%%%%%%%%%%%%%%%%%%
%  Eingabe:       %
%%%%%%%%%%%%%%%%%%%
a = 1;
w = 2*pi;
f = 100;
n = 0.001;
time = (0:n:1)';
time1 = (n:n:1)';
abtast = floor(1/n);
lt = length(time);


%%%%%%%%%%%%%%%%%%%
%  Schwingung     %
%%%%%%%%%%%%%%%%%%%
v_sin = a.*sin(w.*f.*time);
subplot(2,1,1)
plot(time, v_sin)
title('Schwingung')
xlabel('Zeit [s]')
ylabel('Amplitude')
grid on;


%%%%%%%%%%%%%%%%%%%
%  Fourier        %
%%%%%%%%%%%%%%%%%%%
z = fft(v_sin,abtast)
Leistungsspektrum = z.*conj(z)/abtast; %Leistungsdichtespektrum
b = length(Leistungsspektrum);
x = (0:((abtast-1)));
subplot(2,1,2)
stem(x,Leistungsspektrum(1:abtast));
title('Spektrum')
xlabel('Frequenz [Hz]')
ylabel('Amplitude')
grid on;

%%%%%%%%%%%%%%%%%%%%
%  Ausgabe         %
%%%%%%%%%%%%%%%%%%%%
[MaxLeistungsspektrum,freq] = max(Leistungsspektrum);
freqMax = (freq-1);
MaxAmplitude = MaxLeistungsspektrum
MaxFrquenz = freqMax

Autor: die ??? (Gast)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
Doppeltes posting ist unhöflich.

  Beitrag "Fourier Matlab"

Autor: Dennis (Gast)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
Tschuldigung,
letztes Mal kam ziemlich schnell eine Antwort. Da dachte ich, dass die 
vorherige Mail nicht angekommen ist.
Wird nicht wieder vorkommen

Autor: charmin626 (Gast)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
Hallo zusammen

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

Das hat mir sehr weitergeholfen, Matlab zeigt mir jetzt an was ich sehen 
wollte, aber ich würde das ganz gerne auch verstehen. Das meiste ist ja 
bereits erklärt worden, aber was ich nicht verstehe ist wo die 1000 bei 
f=1000*(1:(points/2))/points; herkommt. Wenn mir das jemand kurz 
erklären könnte wäre das super!

Danke!

Autor: Michael Lenz (hochbett)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
Hallo,

> %FFT
> Y = fft(x,points);
> Pyy=Y.*conj(Y)/points; %Leistungsdichtespektrum
> f=1000*(1:(points/2))/points;
> subplot(3,1,2)
> plot(f,Pyy(1:(points/2)));
> grid on;
Zu Deiner Frage: Die 1000 entsprechen der Abtastfrequenz von 1000 Hz, 
die Du mit dem Befehl t = 0:.001:.239 eingestellt hast.

Allerdings enthält der Frequenzvektor noch einen Fehler. Der 
Frequenzvektor muß mit 0Hz anfangen. Er geht dann in Schritten von fs/N 
(fs: Abtastfrequenz, N: Anzahl der Punkte des Zeitvektors) weiter. Bei 
reellwertigen Signalen sind die Spektrallinien mit mehr als fs/2 
redundant.

Richtig ist:
fs = 1000;              % Abtastfrequenz
t = 0:1/fs:.230;        % Zeitvektor
N = length(t);          % Anzahl der Abtastpunkte
f = 0:fs/N:(N-1)/N*fs;  % Frequenzachse

Wenn Du Dir nur die erste Hälfte des Spektrums ansehen willst, mußt Du 
zwischen geradem und ungeradem N unterscheiden. Schau Dir gelegentlich 
auch einmal den Befehl "fftshift" an.


Gruß,
  Michael

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.