Verfasst am: 11.07.2012, 10:39
Titel: projekt chrisl - frequenzanalyse, etc. von scannerpunktwolke
hi leute,
ich mach denn hier mal ein neues thema auf, was ich erstmal etwas allgemeiner formuliere und in dem ich denn bleiben kann^^
ich fummel grad nen code zusammen, der mir ne frequenzanalyse macht. dabei kommen fragen auf - warum muss man das mit dem nullen auffüllen machen? bzw habe ich unterschiedliche ergebnisse, wenn ichs selbst mache oder eben nicht. nachdem ich das mache, sieht die amplitude nach einer normalen fft so unrealistig klein aus.. das ganze ist wohl hauptsächlig ein skalierungsproblem bzw verständnisproblem, wo dabei eine skalierung fehlt bzw falsch ist.. das mit dem fenstern scheint generell nicht all zu gut zu klappen - und bei datensatz x fensterfunction ist das alles bis 2 skaliert. sollte ich für den ausschnitt mit 2 dividieren und anschließend wieder multiplizieren wenn nur das positive sprectrum betrachtet wird?
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)|')
hauptproblematik sind die verschiedenen fft ergebnisse wenn ich nullen selber auffülle oder eben nicht.
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?? ich teste da grad nur
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 Smile 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 (ev. einmal Mi(:,2) durch Yi ersetzen ) 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
mist jetzt habe ich vergessen die datei ranzuhängen. nachtrag:
edit: also ja fensteramplitude kann ich ja selbst skalieren mit 2, allerdings sieht die gefensterte fft komisch aus - warum schluckt er die große amplitude bei 0,2 hz? dürfte eigentlich nicht sein - ist doch kein hochpassfilter.. die frequenz sieht man ja schon analog im rohsignal^^
weiterhin glaube ich, dass mir da n fehler beim nullen auffüllen passiert ist..
Schau dir dort die Skalierung für den Betrag/Amplitude an. Ich empfehle dir zunächst ein einfaches Testsignal (dc + ac*sin(2*pi*f)) zu verwenden, bei dem kein Leakage auftritt. Hier müssen im Betragsspektrum ja auch die richtigen Amplitudenwerte für dc und ac angezeigt werden.
Ist das geschafft, ein Testsignal mit Leakage (f kein ganzes Vielfaches zu df!!) und nun das Fenster einbinden. Hier müsste ja dann eine Besserung gegnüber dem Betragsspektrum ohne Fensterung zu sehen sein. Aber das ist alles in dem Skript in der Skriptecke drin. Warum also alles selber machen, wenn es dort schon ist?
ok danke. ja ich bleibe denn lieber hier in dem thread - ich werde durch antworten nur immer wieder woanders hingelockt^^
guck mir das mal an
und ja ich gucke auch immer alle mögliche skripte an - leider hat mir das bisher noch nie wirklich geholfen - schwierig das auf meins zu übertragen. dagegen, wenn jmd direkt was zu meinen fragen sagt, bringt mich das gleich extrem weiter, daher veruche ich es grad so (was nicht bedeutet, dass ich mir nicht weiter die skriptecke angucke)
gruß
edit: hab mir das beispielskript zum x-mal angeschaut. es ist für einen alten hasen vielleicht kein problem, schnell mal 20 variablen und dimensionen zu ändern um sowas zu verbinden. ich komme so leider nicht weiter. habe schon versucht da was zu integrieren, habe denn aber gemerkt, dass ich deutlich besser vorran kommen und vorallem mehr dabei lerne, wenn ich versuche meinen code weiterzuentwickeln und höchstens mal einzelteile auszuprobiere. so sind mir auch schon einige lichter aufgegangen^^ ich glaube eh, dass die meisten, bevor sie ein thema starten, bereits in die skriptecke geschaut hatten und dabei nicht wirklich weiterkamen.
achja und das mit leakage verfälscht vllt leicht die daten, aber ich glaube nicht, dass das bei meinem beispiel mit 6800 messwerten so ins gewicht fallen kann, aber vllt hab ichs auch nur falsch verstanden.. der 2. schritt vor dem ersten tut mir eh nicht sonderlich gut - ich muss nu gucken wie ich den skalierungsfehler wegbekomme.. und wenn ich eure variante mit meiner vergleiche, sind die jeweils an nem anderen ende falsch, da komme ich nicht so ohne weiteres weiter :S
Zuletzt bearbeitet von Chrislap am 11.07.2012, 12:29, insgesamt einmal bearbeitet
Was du vorher mit deinem Signal machst, ist doch irrelvant. Du kannst für die FFT dir auch eine eigne Funktion schreiben.
Code:
function[mag, mag_dB, fv] = FFT_betragsspektrum( signal, nfft, fa) ;
% Input: % Signal im Zeitbereich % nfft = Anzahl Messwerte für fft % wenn nfft > length(sig) -> fft(sig,nfft) führt Zeropadding durch % fa = Abtastfreq. % Output: % Magnitude des Spektrums linear und dB skaliert % Frequenzvektor fv in [Hz] von 0...fa/2
% un-,gerade Anzahl Messwerte? ifmod(nfft,2) == 0;
k = (nfft/2) + 1;
else
nfft = nfft + 1;
k = (nfft/2) + 1;
end
Als Rückgabewerte erhältst du nun die gewünschten Ergebnisse. Nun kann du mit plot(fv,mag) die lineare Darstellung oder mit plot(fv,mag_dB) die logarithmische (nur Y-Achse) Darstellung machen.
ja aaalso, nu hab ich deine fft-variante genommen - habs aber direkt reingetan, ohne function - das macht es für mich einfacher erstma die ganzen debugs zu finden, die er ausgibt, wenn ich von hier n code reinnehme^^
so sieht das denn gekürzt bei mir aus (nur fft):
Code:
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 (mean) diesen für die Stellen (stelle), 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
L = numel(Mi(:,1)); % Anzahl Werte zur Kontrolle ob Interpolation und Löschen geklappt hat
poti = 2^nextpow2(L); % nächste 2er Potenz / gewünschte FFT-Länge
fa= 50; % Messfrequenz (Eingabewert) ifmod(poti,2) == 0;
k = (poti/2) + 1;
else
nfft = poti + 1;
k = (poti/2) + 1;
end
fn = fa/2; % Nyquistfreq.
df = fa/poti; % Frequenzauflösung des Spektrums % Frequenzvektor: Darstellung bis Nyquistfreq.
fv = 0: df : fn;
sig = Mi(:,2);
% Signal transformieren
Fy = fft(sig,poti);
% Betrag - nur positives Freq.spektrum
Fy_pos0 = abs(Fy(1:k));
% Skalierung
mag = [Fy_pos0(1)/poti; Fy_pos0(2:k-1)/(poti/2); Fy_pos0(k)/poti];
% dB
mag_dB = 20*log10(mag + eps); % eps = kleine Konstante zur Vermeidung von log(0)
fig = figure(1);
plot(fv,mag_dB)
fig = figure(fig+1);
plot(fv,mag)
allerdings gibt es da exakt das gleiche ergebnis wie bei meiner variante - nur das ich bei meiner variante etwas mehr durchblicke^^
ich war da ja schonmal n schritt weiter mit meiner fft und der anderen skalierung (ob durch poti oder L scheint weiterhin egal^^). naja jedenfalls weiß ich jetzt, dass das mit der fft anscheinend gar nix zu tun hat.. so langsam blicke ich nicht mehr durch bei den ganzen varianten wenn ich hier und da mal n anderen code einschiebe..
und die logarytmische darstellung verstehe ich auch nicht so recht - warum ist die immer so im negativen und verzerrt? ist das normal?^^
so also ich habs jetzt wieder mit dem ansatz von dir und heut morgen versucht^^
da bekomme ich anscheinend die richtigen frequenzen für die db_ansicht - vllt kann mir jmd die amplituden erklären^^
allerdings hakt es bei der anderen ansicht - er meint:
Error using plot
Vectors must be the same lengths.
also mg und signal.
Code:
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 (mean) diesen für die Stellen (stelle), 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
Mi(:,1) = Mi(:,1)/50;
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 / gewünschte FFT-Länge
fa= 50; % Messfrequenz (Eingabewert)
a= 1;
%a = poti/L; % ??Amplitudenskalierung??Länge Potenz durch Länge Signal?test
Yi = a*sin(2*pi*fa*Yi);
Yi = [Yi(1:poti/2)zeros(1, poti/2)]; % neuer y_Vektor mit Nullen aufgefüllt
% 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('Single-Sided Amplitude Spectrum of y(t)') % xlabel('Frequency (Hz)') % ylabel('|Y(f)| in cm')
ifmod(poti,2) == 0;
k = (poti/2) + 1;
else
nfft = poti + 1;
k = (poti/2) + 1;
end
amplH = abs(Y); % Betrag bilden -> komplexer Frequenzvektor H für Amplitudengang
df= fa/L;
x_fn = 0 : df : fa/2;
amplitudengang = [amplH(1)/L; amplH(2:L/2)/(L/2); amplH((L/2)+1)/L]; % DC-Bin und Nyquistfreq. auf N normieren! % Plot Ausgabe der Magnitude in dB
plot(x_fn, 20*log10(amplitudengang + eps))% eps = kleine Konstante um log(0) zu vermeiden title('Single-Sided Amplitude Spectrum of y(t)') xlabel('Frequency (Hz)') ylabel('|Y(f)| in dB')
% falsch: hier darf nicht durch poti geteilt werden. Die Skalierung erfolgt lediglich in der Zeile amplitudengang = ...
Y = fft(Mi(:,2),poti)/poti; % FFT
% Um nicht alles umzuschreiben, änder ich hier L % sonst stimmt die Skalierung wieder nicht % und x_fn hat nicht die Länge von amplitudengang, weshalb die Fehlermedlung beim Plot kommt.
L = poti;
df= fa/L;
x_fn = 0 : df : fa/2;
amplitudengang = [amplH(1)/L; amplH(2:L/2)/(L/2); amplH((L/2)+1)/L]; % DC-Bin und Nyquistfreq. auf N normieren!
Ich empfehle dir nach wie vor erstmal mit einem einfachen Testsignal zu arbeiten. Da weißt du was für Amplituden rauskommen müssen und nicht wie bei deinem eigl. Signal. So ist das immer nur Raten, ob das Spektrum richtig ist. Dir scheint es offensichtlich an Grundlagen zu fehlen, siehe dB Darstellung etc. Fang also erstmal einfach an, bevor du hier mit Messdaten arbeitest und versuche das nachzuvollziehen.
ja in der tat kann ich nicht viel mit db-darstellung anfangen - daher hatte ich gefragt - habs jetzt auch auf die schnelle nirgends gefunden. mir ist auch noch nicht bewusst, dass ich das brauche - das kam erst durch euch ins spiel - alles was ich brauche kenne ich ja soweit (jedenfalls die sachen wo ich weiß dass ich sie brauch)^^ ich versuche ja schon immer echt speziell zu fragen, bei den sachen die ich mir nicht selber iwo herholen konnte.
und ja und das mit den durch poti teilen, hatte ich von dir - wohl falsch verstanden und falsch angebracht, aber das macht ja eh keinen merklichen unterschied, hatte ich gesehen^^ deinen ratschlag mit einem bekannten signal werde ich beherzigen - mein signal in einen anderen code zu pressen erscheint mir aufwendig und ich weiß immernoch nicht ob es stimmt^^ außerdem beschleicht mich langsam das gefühl, dass meine auswertung von anfang an gar nicht falsch war - ich habe nämlich mal das ganze so hingebogen, dass ich die ergebnisse bekomme, wie in dem einen code von euch, aber richtig sieht der code jetzt nciht wirklich aus, aber ergebnisse scheinen halt nicht so falsch..
ifmod(poti,2) == 0;
k = (poti/2) + 1;
else
nfft = poti + 1;
k = (poti/2) + 1;
end
amplH = abs(Y); % Betrag bilden -> komplexer Frequenzvektor H für Amplitudengang
df= fa/L;
x_fn = 0 : df : fa/2;
amplitudengang = [amplH(1)/L; amplH(2:L/2)/(L/2); amplH((L/2)+1)/L]; % DC-Bin und Nyquistfreq. auf N normieren!
plot(x_fn, 20*log10(amplitudengang + eps))% eps = kleine Konstante um log(0) zu vermeiden% Plot Ausgabe der Magnitude in dB title('Single-Sided Amplitude Spectrum of y(t)') xlabel('Frequency (Hz)') ylabel('|Y(f)| in dB')
aber so ist dass, wenn ich versuche fremde codes einzubinden und meins verwerfe. da geht mir leider einiges an verständnis verloren.. ich sollte vermutlich den heutigen tag zurückspuhlen und wieder an meinem code versuchen rauszufinden, was da genau das problem beim zeropadding bzw fenstern darstellt..
andere sache: tut mir leid, dass ich hier immer wieder quelltexte reinkopiere - leider kann ich die alten nicht mehr löschen - kann ruhig jmd machen, der die rechte dazu hat..
Aber in diesem Quelltext sind doch schon wieder die gleichen Fehler, die ich in meinem vorherigen Post schon kommentiert habe.
Es macht sehr wohl einen Unterschied, ob man Y=fft()/L bzw. poti teilt, oder das weg lasst. Bei einer Amplitude von 1 kommt dann nämlich nicht die 1 raus sondern 1/L bzw. 1/poti. Das ist doch wohl ein Unterschied.
Ehrlich gesagt verstehe ich wiederrum nicht, was du dort alles mit deinem Messsignal machst. Weder das ständige Interpolieren in der for Schleife noch das Löschen des halben Signals um Nullen anzuhängen, ergibt für mich einen Sinn. Aber nochmal...für die richtige Skalierung ist das erstmal nicht wichtig. Die überprüft man an Daten, wo ich die Ergebnisse bereits kenne.
Die dB Skalierung wird u.a. dazu verwendet, um größere Wertebereich besser darstellen zu können. Ebenfalls ergeben sich andere Verläufe...aus einem exponentiellen Verlauf im Linearen wird eine Gerade durch die Logarithmierung, was für die Analyse von Vorteil sein kann. Es kommt immer drauf an, was man aus dem Spektrum lesen möchte...evtl. reicht ja bei dir die lineare Darstellung.
Edit: Dir sollte ebenfalls bewusst sein, dass durch Interpolieren, Aliasfrequenzen in deinem Signal entstehen, weshalb dann normalerweise eine Aliasingfilter verwendet wird.
ok das gucke ich mir nochmal genauer an. Das ist wohl ein verständnisproblem - die amplitude sollte sich doch aus den daten ergeben und ich dachte es gibt eine übliche skalierung dazu. also du meinst, das teilen sollte man lieber weglassen?! gucke ich mir nochmal an
das interpolieren und löschen hat mit meinen datensatz zu tun, der nicht besonders gut ist. es fehlen ständig werte und denn habe ich doppelte mit dem selben zeitstempel, die leicht abweichen. das hat anscheinend auch ganz gut geklappt. was du mit dem löschen meinst, um nullen ranzuhängen, weiß ich jetzt nicht ganz - das sollte doch keinen großen unterschied machen, da es unter meiner messgenauigkeit liegt - das war auch eher, um bei der rumprobiererei ein wenig zeit zu sparen.. aber das mal beiseite, das ist ja nix großes..
ich hatte halt festgestellt, dass bei der db-variante der frequenzbereich nach rechts gerutscht war und wollte wissen woran das liegt - da muss doch bei einem der beiden varianten ein fehler gewesen sein, der nicht unbedingt ne skalierung ist (die variable amplitudengang ist da wohl der grund)..
ich denke es reicht die lineare darstellung bei meinen miesen daten^^ der letzte punkt ist interessant. es ist klar, dass dadurch alliasfrequenzen auftauchen. allerdings hätte ich irgendwie nicht gedacht, dass man die so gut filtern kann, da die lücken doch so willkürlich auftraten und dann interpoliert wurden. das muss ich mir dann mal genauer angucken.
danke dass du so viel geduld mit mir mitbringst
achja, hast du vllt noch einen vorschlag bzw. gibs hier irgendwo einen geeigneten datensatz, den ich zum probieren nehmen kann? am besten wäre ja auch ne txt mit einigen mehr werten und mehr als eine vorkommende frequenz? danke
fa = 2000; % Abtastfrequenz
fn = fa/2; % Nyquistfrequenz
N = 4096; % 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; % ganzzahliges Vielfaches zu df -> wichtig wg. leakage
f2 = 20.5*df; % hier wird nun leakage auftreten
ac = 1; % Amplitudenvorgabe
dc = 10; % Gleichsignalanteil
y = dc + ac*sin(2*pi*f1*t) + ac*sin(2*pi*f2*t); % y-Vektor
Mehr braucht es erstmal nicht. Somit ist klar, es muss bei 0 Hz im Spektrum eine Ampliutde von 10 dargestellt sein. Bei der Freq. f1 muss die Ampl. = 1 sein. Bei f2 müsste sie theoretisch auch 1 sein, da aber hier leakage autritt, ist so ohne Fenster geringer. Ebenfalls treten um f2 Nebenmaxima auf. Hier kannst du dann noch deine Fensterung testen...da habe ich aber keinen Fehler sehen können. Das sollte also gehen.
Wenn du nun noch Zeropadding testen willst, lass das Anhängen von Nullen doch die fft() Funktion machen.
In deinem Fall wird gilt ja
Code:
poti = 2^nextpow2(L); % L ist die Länge des signals
Was du mit der Verschiebung bei dem anderen Bsp. meinst zwischen dB und lineare Skalierung, kann ich nicht nachvollziehen. Sonst mach doch bitte einen Screenshot.
Ich habe dir hier jetzt mal noch ein kleines Testprogramm geschrieben, an Hand du die Skalierung, Fensterung und Zeropadding überprüfen kannst.
Du musst dazu nur noch meine FFT Funktion von weiter oben in einen seperaten m-file kopieren und mit dem Namen der Funktion, also FFT_betragsspektrum abspeichern.
Vergleiche die Amplitudenwerte des Testsignals (ac und dc mit den Amplituden bei 0 Hz und f1 Hz und f2 Hz) mit den Werten im Spektrum.
fa = 1000; % 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 = 20.5*df; % hier wird leakage auftreten
f2 = 200*df; % ganzzahliges Vielfaches zu df -> wichtig kein leakage
ac = 8; % Amplitudenvorgabe
dc = 10; % Gleichsignalanteil
y = dc + ac*sin(2*pi*f1*t) + ac*sin(2*pi*f2*t); % y-Vektor
% Wechsel in den Freq.-bereich % Betragsspektrum des Testsignals ohne Fensterung
nfft = N; % 1024 Messwerte [mag, mag_dB, fv] = FFT_betragsspektrum( y, nfft, fa ) ;
figure(2) subplot(2,1,1) plot(fv,mag);
xlabel(['Auflösung: ',num2str(df),' Hz Frequenz in [Hz]']);
ylabel('Magnitude') title('Betragsspektrum - ohne Fensterung');
grid on;
% Signal mit Hann Fenster
win = hann(N)';
y_win = y.*win*N/sum(win); % Fensterung mit Amplitudenkorrektur
[mag, mag_dB, fv] = FFT_betragsspektrum( y_win, nfft, fa ) ;
subplot(2,1,2) plot(fv,mag);
xlabel(['Auflösung: ',num2str(df),' Hz Frequenz in [Hz]']);
ylabel('Magnitude') title('Betragsspektrum - mit Hann-Fensterung');
grid on;
% Test - Zeropadding % Es werden nur 900 Messwerte des Testsignals verwendet % 124 Nullen werden am Ende angehängt, um wieder 1024 Werte zu haben
y_900 = y(1:900);
% da nfft immer noch 1024, hängt die Funktion FFT die Nullen automatisch an % denn length(y_900) < nfft [mag, mag_dB, fv] = FFT_betragsspektrum( y_900, nfft, fa ) ;
figure(3) plot(fv,mag); hold on;
xlabel(['Auflösung: ',num2str(df),' Hz Frequenz in [Hz]']);
ylabel('Magnitude') title('Betragsspektrum - Nullen angehängen');
grid on;
% 124 Nullen werden nun selbst angehängt
y_zerop = [y_900 zeros(1,nfft - length(y_900))];
% fft() hängt nun selbst keine Nullen mehr an da length(y_zerop) == nfft [mag, mag_dB, fv] = FFT_betragsspektrum( y_zerop, nfft, fa ) ;
figure(3)
%subplot(2,1,2) plot(fv,mag,'r--');
xlabel(['Auflösung: ',num2str(df),' Hz Frequenz in [Hz]']);
ylabel('Magnitude') legend('fft() an Nullen angehängt','Nullen selbst angehängt');
grid on;
super, danke. das hilft mir weiter. ich muss es mir noch genauer ansehen, aber augenscheinlich scheint mein code soweit okay zu sein, bis auf das die amplitudenskalieren noch nicht stimmt. fenstern scheint ja zu klappen.
danke und gruß, chrisl
edit: und wenn ich so weiter gucke, scheint die amplitudenskalierung auch mit meinem "frequarea" geklappt zu haben. allerdings nur bei meinen daten - bei dem beispielsignal (was übrigens super ist auch um weitere dinge zu testen) hab ich denn iwie ne kommaverschiebung drin
sogar meine andere frage konnte ich mir jetzt etwas genauer selbst erklären. das mit dem nullen ranhängen hatte manuell bei der fensterfunktion doch nicht geklappt. vor der fensterung kann man es eben nicht machen, sonst stimmen dann die dimensionen nicht mehr und danach gehts auch nicht so einfach, weil die kombo aus hann und signal ja sone üble matrix ist^^
soo jetzt habe ich doch noch einmal 2 winzige fragen^^ würde bestimmt in einen anderen abschnitt des forums passen, aber ich machs mal hier - das weiß bestimmt eh jeder von euch - ist ebstimmt ne richtiggehende beleidigung^^
erstens: wie kann ich n bei sublot von einem einzelnen fenster den darstellungsauschnitt ändern - also ne grenze festlegen. bei nem normalen plot hab ichs schon mit axis hinbekommen - bei subplot leider noch nicht..
zweitens: ich will jetzt mein signal erst ab einer bestimmt zeit (x) betrachten, also den anfang quasi abschneiden. so in etwa:
es gibt bestimmt zig varianten einfach ab x in eine neue matrix zu schreiben, weil mit so einer schleifer rattert er sich ja dumm und dämlich - genauso wenn ich ne schleife zum schreiben ab x in eine neue matrix machen^^
gruß
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.
) 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