fa = 8000; % Abtastfrequenz
fn = fa/2; % Nyquistfrequenz
N = 1024; % 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 = 1; % Amplitudenvorgabe
y = a1*sin(2*pi*f1*t); % y-Vektor
y = [y(1:N/2)zeros(1, N/2)];
% Graphische Darstellung % ---------------------- % max. Amplitude zur Skalierung der graphischen Darstellung feststellen:
max_y = max(abs(y))*1.1;
fig = figure(1);
plot(y) axis([0 N -max_y max_y]) title('Datensatz') ylabel('Amplitude') xlabel('N Stützstellen') grid
% 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('Amplitudengang') ylabel('Amplitude') xlabel(['Auflösung: ',num2str(df),' Hz Frequenz in Hz']) grid
% Ausgabe in dB % ------------------
fig = figure(fig+1);
plot(x_fa-fn, 20*log10(amplitudengang))
%axis([-fn fn -10020*log10(a/2)+3]) axis([-fn fn -1003]) title('Amplitudengang') ylabel('Amplitude in dB') 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('Amplitudengang') 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 -10020*log10(a)+3]) title('Amplitudengang') ylabel('Amplitude in dB') xlabel(['Auflösung: ',num2str(df),' Hz Frequenz in Hz']) grid
% Fensterfunktion % ----------------------
% Anhang an die bereits erfolgte Untersuchung % -------------------------------------------
win = hann(N)';
%y_win = y.*win; % Fensterung ohne Amplitudenkorrektur
y_win = y.*win*N/sum(win); % Fensterung mit Amplitudenkorrektur
max_y = max(abs(y_win))*1.1;
fig = figure(fig+1);
plot(y_win) axis([0 N -max_y max_y]) title('Datensatz nach Fensterung mit Hann-Fenster') ylabel('Amplitude') xlabel('N Stützstellen') grid
% Berechnung der FFT % ------------------
H = fft(y_win, 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);
bezüglich zum Script "FFT Umfassendes Beispiel" habe ich ein paar Fragen.
Soll 'df' für 'delta f' stehen und die Fensterlänge angeben, also die Länge des Ausschnittes? Oder was soll die Frequenzauflösung sein? Es ist ja der Quotient aus der Abtastfrequenz und der Anzahl der Abtastungen.
Des weiteren möchte ich anstatt des selbst initialisierten Signals eine kurze , ca. 5 s lange Audiodatei mittels 'wavread'einlesen. Wie würde sich diesbezüglich df und f1 verändern?
Lieben Gruß und Dank
Tim
gast
Gast
Beiträge: ---
Anmeldedatum: ---
Wohnort: ---
Version: ---
Verfasst am: 03.06.2010, 11:59
Titel:
Hallo,
hab eine kleine Frage zu folgendem Teil:
Code:
H = fft(y_win, 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);
danke vorerst mal dass du dein code veröffentlichst. Ich bin an meiner diplomarbeit dran. Ich möchte gerne dein code vewenden, jedoch möchte ich es mit dat.-file verknüpfen. Das Programm soll mein dat-file aus der messung verwenden. Mein datfile enthält 2 zeilen und gewisse Anzahl spalten (so kann man die Abtastfrequenz und die anzahl Messwerte ermitteln). Mein dat.file enhält den Verlauf der Zeit-Amplitude meines Beschleunigungssensors.
Hat jemand eine Lösung wie ich diesen code schreiben kann. Als Grundlage dazu dient der Code von Nils.
ich habe zwei Fragen zur Realisierung der Fensterung.
1.) Frage:
Was steckt hinter der Amplitudenkorrektur? Oder anders gesagt, warum wird (für ein breites Fenster) die Amplitude mit "2" multipliziert?
Anmerkung: sum(Hann(N-> unendlich)) -> N/2
Code:
%y_win = y.*win; % Fensterung ohne Amplitudenkorrektur
y_win = y.*win*N/sum(win); % Fensterung mit Amplitudenkorrektur
2.) Frage:
Für die Fast-Fourier muss ja wegen des Algorithmusses im vgl. zur DFT mit N=2^n Elementen (nextpow2) gearbeitet werden. Entspricht die Anzahl der Elemente nicht diesem Kriterium werden die fehlenden Elemente durch Nullen aufgefüllt.
Muss bei Verwendung der FFT die Fensterung nun VOR oder NACH der "Erweiterung" erfolgen?
Vielen Dank für eine Antwort,
Germanus
Gunnar
Gast
Beiträge: ---
Anmeldedatum: ---
Wohnort: ---
Version: ---
Verfasst am: 01.04.2011, 11:28
Titel:
Moin,
Fensterung immer vor dem ZeroPadding machen. Du willst ja durch die Fensterung vermeiden, künstliche Sprünge zum Signal hinzuzufügen. Wenn du erste Nullen anhängst, und dann das ganze fensterst, hast du irgendwo mittendrin ja immernoch den Sprung.
Gast
Beiträge: ---
Anmeldedatum: ---
Wohnort: ---
Version: ---
Verfasst am: 22.08.2011, 14:06
Titel: Fehlermeldung bei der Analyse von WAVE-Dateien
Hallo Leute,
ich habe das Problem, dass wenn ich anstatt von der Sinus-Halbwelle, die im diesem Beispiel verwendet wird, eine Wave-Datei analysieren möchte, dann bekomme ich eine Fehlermeldung:
??? Error using ==> horzcat
CAT arguments dimensions are not consistent.
Error in ==> fft_t at 84
amplitudengang = [amplH(1)/N amplH(2:N/2)/(N/2)]; % DC-Bin auf N normieren!
Hat jemand eine Ahnung, wo das Problem ist? Vielleicht die Datenmenge oder Eingangsvektorlänge?
Für schnelle Hilfe oder Lösungsvorschläge währe ich dankbar!
ich hätte zu diesem thema auch die eine oder andere frage. und zwar fummel ich grad nen code zusammen, der mir ne frequenzanalyse macht. dabei kommen fragen auf - warum muss man das mit dem nullen auffüllen machen? nachdem ich das mache, sieht die amplitude nach einer normalen fft so unrealistig klein aus.. die befehle dafür habe ich aus diesem thread (ich hab nur variablen ausgetauscht und versucht mit der skalierung was zu machen) --> ist glaube eher n verständnisproblem.
eingebaut sieht das denn so aus bei mir:
I = importdata('1_2500.txt', '\t');
I(:,3) = []; % Höhe löschen, da sie hier nicht gebraucht wird
B = sortrows(I,[12])% Sortieren der Messergebnisse
B(:,1) = round(B(:,1)*10); % profilnummer umdimensionieren disp(num2str(B,'%7.4f')); % runden für weniger rechenaufwand
X = B(:,1); % Matrix zerlegen in X und Y
Y = B(:,2);
xu = unique(X); % schreibe Werte in xu, die kein 2. Mal auftauchen [~, stelle] = ismember(X, xu)% speichere die Stellen, die in X und xu vorkommen
yu = accumarray(stelle, Y, [], @mean)% bilde einen neuen Y-Wert (yu) für jeden X-Wert und mittel diesen für die Stellen (loc), wo es mehrere Y-Werte gibt
for i=1 : max(xu)% Schleife X-Achse für Interpolation, von 1 bis zum maximal auftretenden Wert
Xi = (1:max(xu)); % zähle durch von 1 bis zum letzten Wert
Yi = interp1(xu,yu,Xi,'linear') ; % lineare Interpolation an den Stellen ohne Y-Wert end
Mi = [Xi; Yi]'; % die neue vollständige Matrix
Mi(:,2) = Mi(:,2) - mean(Mi(:,2)); % um den Mittelwert - nur so funzt fft
Mittelwert_Messwerte = mean(round(Mi(:,2))) %Kontrolle
L = numel(Mi(:,1)); % Anzahl Werte zur Kontrolle ob Interpolation und Löschen geklappt hat
AnzahlMessungen_als_Kontrolle = L;
poti = 2^nextpow2(L); % nächste 2er Potenz für FFT
fa= 50; % Messfrequenz (Eingabewert -> gesichert)
a= 1;
%a = poti/L; % ???? Amplitudenskalierung????? Länge Potenz durch Länge Signal?
Yi = a*sin(2*pi*fa*Yi);
Yi = [Yi(1:poti/2)zeros(1, poti/2)]; % neuer y_Vektor mit Nullen aufgefüllt
y_win = Yi.*win*poti/sum(win); % Fensterung mit Amplitudenkorrektur
max_y = max(abs(y_win))*1.1;
fig = figure(fig+1);
plot(y_win)
max_y = max(abs(Yi))*1.1;
axis([0 poti -22]) title('Datensatz nach Fensterung mit Hann-Fenster') ylabel('Amplitude') xlabel('N Stützstellen') grid
Y = fft(y_win,poti)/L; % FFT mit Fensterung
frequArea = fa/2*linspace(0,1,poti/2); % Skalierung des Ergebnisses
fig=(figure(fig+1)) plot(frequArea,2*abs(Y(1:poti/2)))% Plot in Skalierung title('Spektrum nach Fensterung') xlabel('Frequency (Hz)') ylabel('|Y(f)|')
ich habe bisher die erste fft ohne nullen auffüllen also mit M(:,2) - wenn man an dieser stelle Yi nimmt, dann nimmt er es auch schon bei der ersten fft mit rein.. also folgende stelle nimmt er dann mit rein
Code:
a = poti/L; % ???? Amplitudenskalierung????? Länge Potenz durch Länge Signallänge??
a= 1% umgehe meinen test
Yi = a*sin(2*pi*fa*Yi);
Yi = [Yi(1:poti/2)zeros(1, poti/2)]; % neues Y mit Nullen aufgefüllt
aber auch ohne diese stelle, an der ich grade hänge, kann ich nicht so richtig einschätzen, ob das mit der fft und hann-fensterung geklappt hat - sieht ja nicht ganz so toll aus- ich hatte mir erhofft, dass sich die frequenzen mehr von einnander abzeichnen, da der anfang meines datensatz augenscheinlich nicht die beste qualität besitzt.. vielleicht könnte das ja mal einer von euch erfahreren hasen überfliegen ob das richtig sein könnte ich will dabei die frequenz um 0,22 Hz nicht verlieren..
Figure 1: Rohsignal zum Abschätzen der Amplituden
Figure 2: Frequenzanalyse mittels FFT (ohne Nullen)
Figure 3: Hann - Fenster
Figure 4: Dartensatz x Hann-Fenster
Figure 5: neue FFT nach Hann-Fensterung (mit Nullen)
und figure 4 sieht mir auch komisch aus - weiß nicht so recht warum das so genau bis 2 gehen sollte.. kommt halt aus der einen formel weil *poti/sum(win) = 2 ist - aber warum muss man diese amplitudenskalierung anbringen? ich habe den zusammenhang leider nicht ganz gefunden.. liegt der darin, dass nur die hälfte des signals analysiert wird? mit diesem schritt sieht es zumindest meiner ersten fft ähnlich, wenn da auch noch nix mit nullen auffüllen ist..
ich bräuchte halt ne durchgängige auswertemethode mit ähnlichen ergebnissen - ich kann schlecht einmal mit nullen auffüllen und einmal nicht aber bei der gleichen methode passt das nicht so recht zusammen :S
ich lade die beispieldaten mal hoch - ihr bräuchtet nur noch "run" klicken bin für jede hilfe oder auch anregungen was man noch tun könnte, offen. mein gefühl sagt mir, dass da n fehler bei dem nullenauffüllen ist und die fensterung hat bestimmt auch noch nen kleinen haken
danke schön, chrisl
Das Auffüllen von Nullen muss nicht zwangsläufig erfolgen. Ist aber die Anzahl Messwerte nfft (= poti bei dir) keine 2er Potenz, wird der deutlich langsamere DFT Algorithmus verwendet. Die Fast Fouriertransf. (FFT) kann nur bei 2er Potenz erfolgen.
Ist nfft > length(signal) so wird automatisch beim Aufruf von fft(signal,nfft) das Signal auf die Länge nfft durch Nullen erweitert. Dieser Vorgang wird auch Zeropadding genannt. Man kann dadruch auch die Auflösung des Freq.spektrums erhöhen, da sich ja das df = Abtastfreq./nfft ändert. Allerdings fördert dieser Vorgang auch Leakage, nämlich dann, wenn die Signalfreq. kein ganzes Vielfaches von df ist.
Sollten hier noch weitere Fragen, speziell zu deinem Skript sein, würde ich dich allerdings bitten, einen eigenen Thread zu eröffnen. Denn das hat sonst direkt mit dem Skript nichts mehr zu tun.
ok danke - das werde ich denn wohl machen. also ich weiß ja was zeropadding, leakage etc. ist. fachlich verstehe ich das meiste schon und zur not kann ich anchschlagen. ja genau das matlabspezifische - darauf wollte ich abzielen - ob es vom code her nen unterschied macht.. es ist nämlich so, dass ich nen anderes ergebnis bekomme wenn ich mit nullen auffülle (selber).. bin mir unsicher was ich da falsch gemacht haben könnte..
gruß
Hallo,
ich muss auf dem Thema der Amplitudenkorrektur leider weiterhin rumreiten.
Könnte mir jemand erklären, wie der Faktor [code]%N/sum(win)[/code] zustande kommt, sich also herleitet? Hier nochmal der ganze Quellcode:[code]
%y_win = y.*win; % Fensterung ohne Amplitudenkorrektur
y_win = y.*win*N/sum(win); % Fensterung mit Amplitudenkorrektur
[/code]
Vielen Dank
ping
Pepsi
Gast
Beiträge: ---
Anmeldedatum: ---
Wohnort: ---
Version: ---
Verfasst am: 12.09.2012, 15:28
Titel:
hallo
bin als anfänger schrittweise am nachvollziehen der programmierung. komm auch prinzipiell zurecht nur hab ich ein mehr oder weniger triviales problem, an dem ich hängen bleibe...
hier auszüge aus dem code, die ich an das beispiel aus dm forum angepasst habe:
Code:
% Laden des Files
%-----------------------------------
%----------------------------------
%Zeitbereich
% ----------------------------------
x =nums(:,1); % x-Vektor
y =nums(:,2); % y-vektor
Ta = 6.65e-6; % Abtastrate
fa = 1/Ta; % Abtastfrequenz
fn = fa/2; % Nyquistfrequenz
L = length(x); % Länge des Zeitsignals
N = 2^nextpow2(L); % gewünschte FFT-Länge (N=2^x, sonst wird der DFT-Algorithmus verwendet!)
df = fa/N; % Frequenzauflösung
%---Darstellung Datensatz ---
% max. Amplitude zur Skalierung der graphischen Darstellung feststellen:
max_y = max(abs(y))*1.1;
fig = figure(1); %fig
plot(x,y);
%axis([0 N -max_y max_y]) title('Zeitreihe') ylabel('Spannung') xlabel('Zeit') grid
% 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);
ich lese eine excel-datei ein und habe jetzt das problem, bei der fensterung meiner fft folgender fehler auftritt:
Error using .*
Matrix dimensions must agree.
ich hab schon erkannt dass mir in der workspace folgendes angezeigt wird:
y= (250000 x 1 double)
win= (1 x 262144 double)
ich denke, hier liegt der hund begraben, nur steh ich auf dem schlauch, wie ich das problem ohne viel aufwand behebe...!?
und noch zu einer verständnisfrage:
bei der fensterung wird ja die amplitude des frequenzspektrums mit der amplitude des fensters multipliziert, oder ? warum verwendet ihr in eurem beispiel dann bei der fensterung
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.
ich will dabei die frequenz um 0,22 Hz nicht verlieren..
bin für jede hilfe oder auch anregungen was man noch tun könnte, offen. mein gefühl sagt mir, dass da n fehler bei dem nullenauffüllen ist und die fensterung hat bestimmt auch noch nen kleinen haken
RSS