MATLAB: Anfänger möchte FFT erstellen

Computer: Software & Programmierung Programmierung allgemein
#1

Hi Leute,

ich weiß nich mal ob mir einer von euch weiterhelfen kann, wär darüber aber furchtbar dankbar. Ich fang am besten mal an:

Ich mach gerade ein Seminar, bei welchem ich Federkennwerte bestimmen muss. Dazu muss ich auch die Eigenfrequenz der Federn bestimmen.
Nun hab ich 50k Werte mit einer Abtastrate von 20k Hz ermittelt mit den zugehörigen Zeitangaben. Dabei wird die Kurve aber von einem störenden Rauschen überlagert. Dieses kann ich laut meinem Seminarbetreuer mit Hilfe von matlab und einer FFT eliminieren.

Ich hab jetzt schon einiges dazu durchgelesen, aber als kompletter noob bekomm ich einfach nichts gebacken. Hab auch schon versucht anhand von Beispielen die Programmierung nachzuvollziehen, doch selbst die sind mir zu komplex.

Bis jetzt hab ichs nur geschaftt die Daten mit xlsread einzulesen und zu plotten. Jetzt müsste ich doch schon die FFT anwenden können, oder? Bloß weiß ich nicht wie ich das anstelle.

Könnt ihr mir evtl eine kleine Starthilfe geben, wie ich anfangen könnte oder ob ich überhaupt auf dem richtigen Weg bin?

Entschuldigt bitte meine laienhafte und unpräzisse Ausdrucksweise, ich bin leider nicht vom Fach.

Über jede Hilfe wär ich tausendfach dankbar!

#2

Hallo,

ich würde doch eher dazu raten, Kommilitonen zu fragen.

Es gibt bei Matlab eine Funktioni FFT. Hilfe kriegst Du immer mit z.B.

help fft

oder

doc fft

Da findest Du extrem viel Info. Evtl kannst Du die Schwingung zeit_bereich (die Du mit xlsread gelesen hast) in den Frequenzbereich wandeln mit

f_bereich = fft(zeit_bereich)

Im Frequenzbereich müsstest Du dann deutliche Peaks erkennen. Alle Werte, die keine Peaks sind (also kleiner als threshold), sind Rauschen. Die kannst Du zu Null setzen mit

f_bereich(f_bereich

#3

Hallo,

du kannst über eine FFT die Frequenz schätzen. Hier ist ein kurzes Beispiel:

clear

% Parameter
Nfft = 1024; % Anzahl Abtastwerte für die FFT
fa = 20000; % Abtastfrequenz
fmax = 0.5*fa; % maximale Frequenz die ohne Aliasing aufgelöst werden kann

% Testsignal berechnen
f = 3000; % Sinus-Frequenz
t = (0:frowning:Nfft-1)) / fa; % Zeitachse
x = sin(2*pi*f*t); % Sinus mit Frequenz f
n = randn(size(t)); % Rauschsignale
s = x + n; % Testsignal

% Testsignal plotten
figure(1); clf;
subplot(2,1,1);
plot(t, x);
xlabel(‚Zeit‘);
ylabel(‚Signalwert‘);

% Absolutwert der FFT über der Frequenz plotten
subplot(2,1,2);
sfft = abs(fft(s));
plot(fa * (0:frowning:Nfft/2-1)) / Nfft, sfft(1:Nfft/2));
xlabel(‚Frequenz‘);
ylabel(‚Amplitude‘);

In der FFT siehst du ein Maximum bei der Frequenz des Sinus. Du kannst jetzt das Signal s durch die ersten Nfft Abtastwerte deines Messsignals ersetzen.

Je mehr Abtastwerte deines gesamten Signals du nimmst, desto besser ist die Auflösung der Frequenz, aber desto höher ist auch der Rechenaufwand (für eine einmalige Berechnung mit Matlab natürlich kein Problem).

Die Theorie dahinter ist etwas aufwändiger, wenn du Fragen hast nur zu.

Gruß
Markus

#4

Ahoi,

also wie ich das sehe, ist das sogar eine perfekte Anwendung für das Matlab-Beispiel aus der Hilfe (Befehl: doc fft). Ich nehme zwar an, dass der eigentliche Versuch etwas komplexer war, doch stelle ich mir mal ein einfaches Federpendel vor.

clear all;
% Die folgend Angaben wurden aus der fft-Hilfe entnommen und angepasst
Fs = 20000; % Sampling frequency -> Deine Abtastrate von 20 kHz
T = 1/Fs; % Sample time
L = 50000; % Length of signal -> Die Anzahl deiner abgetasteten Werte

%Den Zeitvektor und den Wertevektor hast du ja schon durch xlsread
%erhalten. Übergib also einfach t=Deine Zeitwerte und y=Deine Messwerte.
%Ich baue mir mal ein Signal zusammen

f_Feder=50; %Ja! Meine Feder hat eine Eigenfrequenz von 50 Hz

t = (0:L-1)*T; % Time vector -> 50k Werte im Abstand von T
% Hier bastel ich mir mal einfach ein Signal, wie ich es von einer Feder
% erwarte.
x = 0.010*sin(2*pi*f_Feder*t); % Auslenkung: 1 cm
y = x + 0.05*randn(size(t)); % Sinusoids plus noise
plot(t(1:50000),y(1:50000)) % Ich hoffe, dein Zeitsignal sieht ähnlich aus, sonst lohnt sich die FFT nicht ^^
figure(1);
title(‚Zeitsignal‘)
xlabel(‚Zeit [s]‘)
ylabel(‚Federposition [m]‘)

NFFT = 2^nextpow2(L); % Next power of 2 from length of y
Y = fft(y,NFFT)/L;
f = Fs/2*linspace(0,1,NFFT/2+1);

% Plot single-sided amplitude spectrum.
figure(2)
plot(f,2*abs(Y(1:NFFT/2+1))) %Hier zeigt er dir nun das ganze Spektrum an
%Nun guckst du dir an, wo das Maximum ist durch heranzoomen oder so
title(‚Single-Sided Amplitude Spectrum of y(t)‘)
xlabel(‚Frequency (Hz)‘)
ylabel(’|Y(f)|’)
%Oder du lässt es dir von Matlab ausgebn
[a_max, i_max]=max(Y(1:NFFT/2+1));
disp([‚Bei einer Frequenz von ‚,f(i_max),‘ Hz liegt die Eigenfrequenz der Feder‘]);

#5

Vielen Dank, du hast mir unheimlich weitergeholfen!
Habs mit deiner Hilfe hinbekommen, wirklich sehr gut erklärt.

Ich wünsch dir ne schöne Woche!

#6

Ich bins nochmal. Ich dachte nur, es hätte geklappt. Ich hab nämlich keine wirklichen Peaks herausbekommen. Nur einen in der nähe von Null. Hier hab ich nen Screenshot: http://img714.imageshack.us/i/frequenzunskaliertdiag…
und hier die zugehörige Programmierung:
http://img198.imageshack.us/i/frequenzunskaliertprog…
Hast du eine Ahnung, was da schief gelaufen ist?

Grüße und nochmals Vielen Dank

#7

Vielen Dank erstmal. Ich hab versucht nach deiner Anleitung vorzugehen, hab Sie auch verstanden.
Nur hab ich dann in etwa so etwas herausbekommen:
http://imageshack.us/photo/my-images/714/frequenzuns…

ich hab also keine erkennbaren Peaks bekommen. Kann das sein, oder hab ich etwas falsch gemacht?

Grüße und nochmals vielen Dank!

#8

Ahoi,

tatsächlich kann man hier recht wenig erkennen. Das liegt jedoch daran, dass das Bauteil nicht von Anfang an schwingt, sondern dass die Erregung noch in den Messwerten aufgenommen wurde. Mit anderen Worten: Zu Anfang ist das Teil in Ruhe - also kein Signal ist vorhanden. Im Zeitbereich von 0 bis etwa 0,3 s hast du also nur einen Gleichanteil von 6. Danach wird es von etwa Zeitindex 0,3 bis 0,6 s stark ausgelenkt und losgelassen (das wäre eine Halbschwingung). Dies geht ganz stark in die FFT ein. Ohne dabei sicher zu sein, behaupte ich nun einfach mal, dass du allein durch das Auslenken eine hohe Amplitude in der Frequenz von 1/(2*0,3 s), also etwa 2 Hz bekommst, sowie weitere Amplituden bei Vielfachen von dieser Frequenz.

Tipp 1: Schneide die Erregung des Bauteils ab. Nimm also nur das Signal ab dem ersten „Nulldurchgang“ (Du hast einen Gleichanteil von 6, also wäre „Sechsdurchgang“ das bessere Wort) bei etwa 0,6 s. Dann hast du nur noch die Schwingung als Signal. Bei 50 000 Messwerten dürfte dies etwa ab dem 12 000. Wert der Fall sein. Das kannst du ganz leicht überprüfen mit plot(t(12000:50000),y(12000:50000); Anschließend legst du dir einen Wertevektor mit den dargestellten Werten an, und machst damit die Fouriertransformation. Also yhilf=y(12000:50000); clear y; y=yhilf; Dadurch verändert sich auch die Anzahl deiner abgetasteten Werte L=50000-12000; Wenn ich nichts vergessen habe solltest du dann weitermachen können, wie beschrieben.

Tipp 2 (optional): Es bietet sich an, den Gleichanteil herauszunehmen, da du nur die Eigenfrequenz haben möchtest. Dies ginge ganz einfach mit y=y-6; vor der Transformation. Es ist nicht wichtig. Nur so vermeidest du den Ausschlag bis zu einer Amplitude von 30000 bei der Frequenz 0 (und andere Peaks fallen ggf. besser auf, sodass du nicht die Lupenfunktion nutzen musst). (Der Wert für die Amplitude müsste eigentlich bei 12 liegen, aber die Amplitude ist vorerst nicht wichtig, da uns ja nur die Frequenz interessiert. Mein Tipp dafür ist, dass die Amplituden nicht durch die Länge des Signals geteilt wurden)

Tipp 3: Zoom. Besonders wenn du noch den Gleichanteil drin hast, muss du das Signal genauer unter die Lupe nehmen.

Damit sollte das dann hinhauen.

Gruß, Ulfson

#9

Nachtrag: Habe mal bisschen gebastelt ^^

clear all;

% Die folgend Angaben wurden aus der fft-Hilfe entnommen und angepasst
Fs = 20000; % Sampling frequency -> Deine Abtastrate von 20 kHz
T = 1/Fs; % Sample time
L = 50000; % Length of signal -> Die Anzahl deiner abgetasteten Werte

%Den Zeitvektor und den Wertevektor hast du ja schon durch xlsread
%erhalten. Übergib also einfach t=Deine Zeitwerte und y=Deine Messwerte.
%Ich baue mir mal ein Signal zusammen

f_Feder=1000; %Ja! Meine Feder hat eine Eigenfrequenz von 1000 Hz

t = (0:L-1)*T; % Time vector -> 50k Werte im Abstand von T
%% Hier bastel ich mir mal einfach ein Signal, wie ich es eben gesehen habe.
x(1:6000) = 0.06;
x(6001:11999)=0.06-2/3*(0.00005.*(1:5999)).^2;
x(12000:50000)=0.06+exp(-(1.4*t(12000:50000)).^2).^5.*sin(2*pi*f_Feder*t(12000:50000));
y = x + 0.0005*randn(size(t)); % Sinusoids plus noise
figure(1);
plot(t(1:50000),y(1:50000)) % Ich hoffe, dein Zeitsignal sieht ähnlich aus, sonst lohnt sich die FFT nicht ^^

title(‚Zeitsignal‘)
xlabel(‚Zeit [s]‘)
ylabel(‚Federposition [m]‘)

%% FFT für alles

NFFT = 2^nextpow2(L); % Next power of 2 from length of y
Y = fft(y,NFFT)/L;
f = Fs/2*linspace(0,1,NFFT/2+1);

% Plot single-sided amplitude spectrum.
figure(2)
plot(f,2*abs(Y(1:NFFT/2+1))) %Hier zeigt er dir nun das ganze Spektrum an
%Nun guckst du dir an, wo das Maximum ist durch heranzoomen oder so
title(‚Single-Sided Amplitude Spectrum of y(t)‘)
xlabel(‚Frequency (Hz)‘)
ylabel(’|Y(f)|’)

%% Ausblendung der unwichtigen Signale und dann FFT
L=50000-12000;
yhilf=y(12000:50000);
clear y;
y=yhilf-0.06; %durch -0.06 wird der Gleichanteil eleminiert danach ist alles gleich

NFFT = 2^nextpow2(L); % Next power of 2 from length of y
Y = fft(y,NFFT)/L;
f = Fs/2*linspace(0,1,NFFT/2+1);

% Plot single-sided amplitude spectrum.
figure(3)
plot(f,2*abs(Y(1:NFFT/2+1))) %Hier zeigt er dir nun das ganze Spektrum an
%Nun guckst du dir an, wo das Maximum ist durch heranzoomen oder so
title(‚Single-Sided Amplitude Spectrum of y(t)‘)
xlabel(‚Frequency (Hz)‘)
ylabel(’|Y(f)|’)
%Oder du lässt es dir von Matlab ausgebn
[a_max, i_max]=max(Y(1:NFFT/2+1));
disp([‚Bei einer Frequenz von ‚,num2str(f(i_max)),‘ Hz liegt die Eigenfrequenz der Feder‘]);

#10

Hallo,

der Peak bei Null liegt daran dass dein Signal nicht mittelwertfrei ist, du solltest diesen also zuerst abziehen:

x2 = x - mean(x);

Außerdem solltest du nur den eingeschwungenen Teil des Signals ab ca. 0,7 Sekunden auswerten. Das Einschwingen ist nicht sinnvoll durch eine Frequenz zu beschreiben.

Du kannst ja dein Ergebnis noch mal hochladen.

Gruß
Markus

#11

Hallo,

eine FFT ist einfach eine Fouriertransformation. Damit kannst Du Schwingungsfrequenzen ermitteln. Du brauchst dazu den Ort des schwingenden Objekt als Funktion der Zeit, also genau das was Du als Input-Daten hast.
Angenommen dieses array von Werten hat den Namen X. Du kannst diese Datenreihe dann in Matlab folgendermassen transformieren:
Y = fft(X)
anschliessend kannst Du Dir das Ergebnis mit plot(abs(Y)) anschauen. Du wirst vermutlich ein Rauschen mit einer deutlich herausragenden Spitze sehen. Die Position dieses Peaks gibt Dir die Eigenfrequenz an.

Wenn Du es genauer wissen willst müsste ich ein bischen ausholen. Dafür brauchst Du aber einige mathematische Kenntnisse. Ggf. einfach nochmal melden.

Viele Grüße,

Christoph

#12

Danke Christoph! Wenn ich noch Fragen habe, melde ich mich.
Viele Grüße,
Stefan

#13

Sau cool, danke. Hab die Antwort von dir erst heute bemerkt, hab das irgendwie übersehen. Probier ich am Dienstag gleich aus und geb dir Bescheid.
Ihr seid mir eine riesige Hilfe, vielen Dank!