könnte sich vielleicht jemand meinen Code anschauen, würde den gerne so effektiv wie möglich haben, damit er in möglichst kurzer Zeit abgearbeitet werden kann!!
Hauptsächlich gehts mir eigentlich um die for - Schleife, wenn sonst noch jemandem was auffällt bitte sagen!!
Im Prinzip werden zwei vektoren komponentenweise verglichen und immer der größte wert in einen neuen Vektor übernommen!!
Danke schon mal für eure Mühen
Code:
%% Vektoren Vergleichstest
% vek1 = 10.*randn(1,10); % vek2 = 10.*randn(1,10); closeall tic
f1 = 10; % signal frequency1
f2 = 20;
fs =200; % sampling frequency
L = 4096; % Signallänge
Ts = 1/fs; % time interval
t = (0:L-1)*Ts; % time vector
y1 = 2*sin(2*pi*f1*t);
y2 = 0.5*sin(2*pi*f2*t);
amax = max(y1)*1.1;
plot(t,y1,'b', t, y2,'r');
axis([01 -amax amax]);
Du hast aber noch einen kleinen Fehler bei der Skalierung der Amplitudenwerte. Du stellst ja nur den positiven Frequenzbereich dar und dann werden der Gleichsignalanteil Y1(1) bzw. Y2(1) und bei der Nyquistfrequenz Y1((L/2) + 1) nicht mit 2 multipliziert.
Hab deine Anmerkung angepasst!
Hab die Vektoren jetzt logarithmisiert, ist das starke auseinandertriften des Spektrums nach den Peaks aufgrund des Leakage zu erklären??
Da meine Perioden kein ganzzahliges Vielfaches meines FFT - Fensters sind??
Hab deine Anmerkung angepasst!
Hab die Vektoren jetzt logarithmisiert, ist das starke auseinandertriften des Spektrums nach den Peaks aufgrund des Leakage zu erklären??
Da meine Perioden kein ganzzahliges Vielfaches meines FFT - Fensters sind??
Das ist definitiv Leakage, da die Frequenzauflösung auch in diesem Bsp. kein ganzzahliges Vielfaches von f1 und f2 ist.
df = fs/L
f1/ df = 1048.6
f2/ df = 1667.7
Was meinst du denn mit gerade und ungerade...die Anzahl an Messwerten? Du solltest schon eine 2er Potenz haben, denn ansonsten wird mit der viel langsameren dft gerechnet.
Die Normierung ist auf jeden Fall richtig. Warum es in dem Mathworks Bsp. falsch steht, kann ich dir nicht beantworten. Es lässt sich ja ganz einfach überprüfen. Wenn du zu y1 z.B. einen Gleichanteil von 10 addierst, muss bei f = 0 Hz ja auch eine Amplitude von 10 stehen. Nach dem Mathworks Bsp. kommen dann aber 20 raus und das ist ja falsch. Ob man die Nyquist- oder Shannonfreq. f((L/2)+1) darstellen sollte, ist ohnehin fraglich. Dieser Wert ist nur selten korrekt berechnet. Aber auch hier darf nicht mit 2 multipliziert werden.
Eine log. Darstellung der x-Achse macht hier aber wenig Sinn. Das ist ja nur dafür gedacht, um große Frequenzbereich besser/übersichtlicher darstellen zu konnen. Bei einem 100 Hz Bereich ist das ja nicht der Fall.
Ok, hatte ich also richtig vermutet, dass das leakage ist!! Danke!
Hab weiter an meinen Programm gearbeitet und eben jetzt auch die Segmentierung eingefügt, um meine geforderte PEAK Darstellung zu realisieren!
Doch obwohl jetzt schrittweise gefenster wird ist der leakage hier fast noch größer als bei allen anderen plots und ich glaub die Amplituden stimmen auch nicht!
(Wann braucht man Amplitudenkorrektur bei Fensterung und wann nicht?
wie hier:
Zitat:
FFT: Umfassendes Beispiel
Code:
%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;
Könntest du dir mein Programm nochmals durchsehen und mir erklären was ich noch besser machen kann und wo mein Fehler liegt???
Code:
%% Vektoren Vergleichstest
% vek1 = 10.*randn(1,10); % vek2 = 10.*randn(1,10); closeall clear tic
f1 = 10.0000; % signal frequency1
f2 = 20.0000; % signal frequency2
fs = 200; % sampling frequency
L = 4096; % Signallänge
Ts = 1/fs; % time interval
t = (0:L-1)*Ts; % time vector
y1 = 2*sin(2*pi*f1*t);
y2 = 0.5*sin(2*pi*f2*t);
amax = max(y1)*1.1;
% graphics
Signal_Plot = figure;
plot(t,y1,'b', t, y2,'r');
grid on;
axis([01 -amax amax]);
%% FFT der beiden Sinus Signale
Y1 = fft(y1)/L;
Y2 = fft(y2)/L;
f = fs/L*(0:L-1); % Frequenzvektor wird gebildet % graphics
FFtplot = figure;
Y1abs = abs(Y1(1:(L/2)+1)); % Absolutwerte der komlpexen FFT - Daten
Y2abs = abs(Y2(1:(L/2)+1));
Y1abs(2:end-1) = 2*Y1abs(2:end-1); % Gleichanteil und Anteil bei Nyquistfrequenz nicht mit 2 multipliziert
Y2abs(2:end-1) = 2*Y2abs(2:end-1);
plot(f(1:(L/2)+1),20*log10(Y1abs),'b', f(1:(L/2)+1), 20*log10(Y2abs),'r');
grid on;
%% Nur die Maximum Werte der jeweiligen Vektoren übernehmen
neuVek = max(Y1abs,Y2abs); % Max - Werte der jeweiligen Vektoren werden übernommen! % Plot des maximalen Vektors
FFtNeuVekPlot = figure;
%% Die zwei Sinussignale werden miteinander addiert
x = y1 + y2;%+2.7653*randn(size(y1)); % Addieren der beiden Sinussignale
Yx = fft(x)/length(x); % FFT der addierten Signale
Add_Plot = figure;
subplot(2,1,1);
Yxabs = abs(Yx(1:(L/2)+1));
Yxabs(2:end-1) = 2*Yxabs(2:end-1); % Gleichanteil und Anteil bei Nyquistfrequenz nicht mit 2 multipliziert
%% Segmentweiser Fensterung mit Overlab und FFT % *************************************************************************
nfft = L/8; % nfft entpricht der Segmentlänge, kleines nfft --> feinere Frequenzauflösung
overlap = nfft/2; % Überlappung zweier aufeinanderfolgender Fenster; nfft = 512
win = hann(nfft,'periodic'); % Hanning - Fenster zur Bearbeitung % Anzahl der zu bearbeitenden Segmente feststellen
k = ((L - overlap)/(length(win) - overlap));
H = zeros(nfft,1); % Ergebnismatrix k. Segmente der Länge nfft transformieren
H_pos0 = zeros(length(0:(nfft/2)),k);
m = 1; % Indexstartwert für das Segment
x = x'; % zu Spaltenvektor transponieren, da win Spaltenvektor ist
ff = 1; % Laufindex für Peak - Finding for ii=1:k
% Fensterung mit Hann window % entnehme aus Signal das Segment der Länge nfft = 512 und multipliziere mit dem Hann-Fenster
signal_win = x(m:nfft+m-1) .* win;
% Segment transformieren
H(1:nfft) = fft(signal_win, nfft)/nfft;
% aus H werden nur die Werte 1...257 = (nfft/2) + 1 entnommen
H_pos0(1:(nfft/2)+1,ii) = abs(H(1:(nfft/2)+1));
ifisequal(ii,2)
spec_peak = max(H_pos0(:,1), H_pos0(:,2)); % max Werte in spec_peak
spec_peak = max(spec_peak,H_pos0(:,ii)); % bei jedem Schleifendurchlauf wird spec_peak upgedated
end % Indexstartwert für nächstes Segement bestimmen
m = m + nfft - overlap;
end
%% Plotten der PEAK Detektion
% Auflösung des Spektrums
f_peak = fs/nfft*(0:nfft/2);
PEAK_PLOT = figure;
plot(f_peak, 20*log10(spec_peak)','b');
title('The PEAK - Plot of a mixed Sine Wave');
xlabel('Frequenz in [Hz]');
ylabel('Amplitude in dB');
grid on;
Ich habe mal noch ein paar Sachen eingefügt. Du solltest dir mal besonders den letzten Plot anschauen und dich dann fragen...macht die Segmentierung und der Detektor Sinn? Aber ich verstehe die Thematik eh nicht so ganz. Auch nicht mit dem neuen Skript, welches du noch gepostet hattest. Ich kenne mich auch nicht in Nachrichtentechnik aus...ich komme aus dem Bereich Automatisierungs- und Informationstechnik
Bei spec_peak hattest du übrigens die Normierung vergessen, weshalb die Amplitude ja nicht stimmen konnte. Ich habe die Skalierung in dB beim letzten Plot auch mal weggelassen, damit man nicht erst umrechnen muss und sie gleich mit den tatsächlichen Amplituden vergleichen kann. Außerdem sollte dir bewusst sein, dass durch Leakage nicht nur das Linienspektrum um die Signalfreq. verbreitert ist, sondern auch einen Fehler in der Amplitude erzeugt. Dieser Fehler ist genau dann maximal, wenn die Signalfrequenz genau in der Mitte von df liegt. Daran ändert auch dann eine höhere Auflösung nichts.
Die Normierung des Fensters solltest du dann nutzen, wenn du an dem Amplitudenwert interessiert bist. Noch was zum Thema Fenster. Soll das Linienspektrum bei Leakage um die Signalfreq. vermindert werden, wählt man z.B. Hann oder Hamming. Hier liegt das Augenmerk dann auf der Frequenz. Wenn man allerdings an einer genaueren Ampl. interessiert ist, wählt man ein Flattopwin. Da wird der Leakageffekt nicht so gut minimiert wie bei Hamming aber die Ampl. ist meist genauer.
Code:
%% Vektoren Vergleichstest
% vek1 = 10.*randn(1,10); % vek2 = 10.*randn(1,10);
%closeall clear tic
f1 = 10.0000; % signal frequency1
f2 = 20.0000; % signal frequency2
fs = 400; % sampling frequency
L = 4096; % Signallänge
Ts = 1/fs; % time interval
t = (0:L-1)*Ts; % time vector % einfacher Test für Leakage % Test = 0 -> kein Leakage
df = fs/L % Frequenzauflösung des Spektrum
test = mod(f1,df)
test = mod(f2,df)
%% FFT der beiden Sinus Signale
Y1 = fft(y1)/L;
Y2 = fft(y2)/L;
f = fs/L*(0:L-1); % Frequenzvektor wird gebildet % graphics
FFtplot = figure;
Y1abs = abs(Y1(1:(L/2)+1)); % Absolutwerte der komlpexen FFT - Daten
Y2abs = abs(Y2(1:(L/2)+1));
Y1abs(2:end-1) = 2*Y1abs(2:end-1); % Gleichanteil und Anteil bei Nyquistfrequenz nicht mit 2 multipliziert
Y2abs(2:end-1) = 2*Y2abs(2:end-1);
plot(f(1:(L/2)+1),20*log10(Y1abs+eps),'b', f(1:(L/2)+1), 20*log10(Y2abs+eps),'r');
grid on;
%% Nur die Maximum Werte der jeweiligen Vektoren übernehmen
neuVek = max(Y1abs,Y2abs); % Max - Werte der jeweiligen Vektoren werden übernommen! % Plot des maximalen Vektors
FFtNeuVekPlot = figure;
%% Die zwei Sinussignale werden miteinander addiert
x = y1 + y2;%+2.7653*randn(size(y1)); % Addieren der beiden Sinussignale % Fensterung für komplettes Signal
win = hann(length(x),'periodic');
Yx = fft(x.* win'*length(x)/sum(win))/length(x); % FFT der addierten Signale
Add_Plot = figure;
subplot(2,1,1);
Yxabs = abs(Yx(1:(L/2)+1));
Yxabs(2:end-1) = 2*Yxabs(2:end-1); % Gleichanteil und Anteil bei Nyquistfrequenz nicht mit 2 multipliziert
plot(t,x);
amaxs = max(x)*1.1;
axis([01 -amaxs amaxs]);
xlabel('Time in [s]');
ylabel('Amplitude');
grid on;
% FFT plot mixed Signal!! subplot(2,1,2);
semilogx(f(1:(L/2)+1), 20*log10(Yxabs+eps),'r'); % eps ist eine kleine Konstante um log10(0) zu verhindern grid on;
%% Segmentweiser Fensterung mit Overlab und FFT % *************************************************************************
nfft = L/8; % nfft entpricht der Segmentlänge, kleines nfft --> feinere Frequenzauflösung
overlap = nfft/2; % Überlappung zweier aufeinanderfolgender Fenster; nfft = 512
win = hann(nfft,'periodic'); % Hanning - Fenster zur Bearbeitung % Anzahl der zu bearbeitenden Segmente feststellen
k = ((L - overlap)/(length(win) - overlap));
H = zeros(nfft,1); % Ergebnismatrix k. Segmente der Länge nfft transformieren
H_pos0 = zeros(length(0:(nfft/2)),k);
m = 1; % Indexstartwert für das Segment
x = x'; % zu Spaltenvektor transponieren, da win Spaltenvektor ist
ff = 1; % Laufindex für Peak - Finding for ii=1:k
% Fensterung mit Hann window % entnehme aus Signal das Segment der Länge nfft = 512 und multipliziere mit dem Hann-Fenster
signal_win = x(m:nfft+m-1) .* win*nfft/sum(win);
% Segment transformieren
H(1:nfft) = fft(signal_win, nfft)/nfft;
% aus H werden nur die Werte 1...257 = (nfft/2) + 1 entnommen
H_pos0(1:(nfft/2)+1,ii) = abs(H(1:(nfft/2)+1));
% Normierung fehlte
H_pos0(2:(nfft/2),ii) = 2 * H_pos0(2:(nfft/2),ii);
ifisequal(ii,2)
spec_peak = max(H_pos0(:,1), H_pos0(:,2)); % max Werte in spec_peak
Das "ii ~= 1" ist überflüssig, weil das für alle Zahlen gilt, wenn sie >= 3 sind.
Die IF-Abfragen innerhalb der Schleife kosten Zeit. Es lohnt sich, solche Ausnahmen einfach explizit vor die Schleife zu ziehen:
Code:
nfft = L/8; % nfft entpricht der Segmentlänge, kleines nfft --> feinere Frequenzauflösung
overlap = nfft/2; % Überlappung zweier aufeinanderfolgender Fenster; nfft = 512
win = hann(nfft,'periodic'); % Hanning - Fenster zur Bearbeitung % Anzahl der zu bearbeitenden Segmente feststellen
k = ((L - overlap)/(length(win) - overlap));
H = zeros(nfft,1); % Ergebnismatrix k. Segmente der Länge nfft transformieren
H_pos0 = zeros(length(0:(nfft/2)),k);
m = 1; % Indexstartwert für das Segment
x = x'; % zu Spaltenvektor transponieren, da win Spaltenvektor ist
ff = 1; % Laufindex für Peak - Finding
c1 = (nfft/2)+1;
ii = 1;
signal_win = x(m:nfft+m-1) .* win;
H = fft(signal_win, nfft)/nfft;
H_pos0(1:c1,ii) = abs(H(1:c1));
m = m + nfft - overlap;
ii = 2;
signal_win = x(m:nfft+m-1) .* win;
H = fft(signal_win, nfft)/nfft;
H_pos0(1:c1,ii) = abs(H(1:c1));
spec_peak = max(H_pos0(:,1), H_pos0(:,2));
m = m + nfft - overlap;
for ii = 3:k
signal_win = x(m:nfft+m-1) .* win;
H = fft(signal_win, nfft)/nfft;
H_pos0(1:c1,ii) = abs(H(1:c1));
spec_peak = max(spec_peak,H_pos0(:,ii));
m = m + nfft - overlap;
end
Ok ich werd mir das jetzt gleich mal alles ansehen und wahrscheinlich(sogar ziemlich sicher werd ich noch ein paar Fragen haben, hoffe des geht in Ordnung )
Ich weiß das die Segmentierung mehr Fragen aufwirft als sie beantwortet, aber es soll halt mal so gemacht werden, werd da auf alle Fälle nochmal nachhaken! Aber mir wurde es halt mal so erklärt!!
@ Jan:
Danke für deine Anmerkung!
Und diese Anordnung ist schneller???????
Aber es ist gut das du es erwähnst, denn mein späteres(End-) Programm sollte so effektiv wie möglich sein, da ich später große Datenmengen zu verarbeiten haben werde!
Allgemein:
Ist das eigentlich in Ordnung, wenn ich immer wieder fertige Programmteile zum " Korrekturlesen" poste, oder verfälscht das den Sinn des Forums???
Ob das schneller ist, weiß ich nicht. Aber Du kannst es einfach mal per TIC/TOC messen.
Es ist im Allgemeinen wichtig, die tatsächliche Größe der Daten zu kennen. Ein Code, der für 5kB Daten extrem effizient ist, weil alle Daten in den prozessor-cache passen, kann bei 5GB Daten gnadenlos scheitern. Im Forum findet man immer wieder Angaben wie "viele Daten" oder "ein sehr großes Array". Das hilft aber nicht weiter, denn Programmieranfänger finden eine 1000*1000 Matrix schon recht groß, für ein FEM-Modell eines PKWs ist das dagegen winzig.
Das Posten von Code und Fragen nach Verbesserungen entspricht dem Sinn des Forums.
Die Fragen zur Fensterung/Frequenzspektrum etc. sind auch beantwortet?! Das ging ja diesmal fix
Einstellungen und Berechtigungen
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.
)
RSS