Hab den Code erstmal 1:1 kopiert und getestet, scheint zu laufen.
Dann habe ich versucht diesen Code auf meinen Signal zu "überstülpen".
Leider scheint dies nicht zu funktionieren, vielleicht kann der eine oder andere mir helfen.
Was Matlab angeht, bin ich ein Anfänger.
Ich muss bei meiner Aufgabe die Grundwelle des gemessenen Signals filtern.
Aus dieser Grundwelle soll dann die Effektivwerte und später dann die Grundschwingungsgehalt sowie Oberschwingungsgehalt berechnet werden.
Bei t= zeit und y = current habe ich meine Vektoren angefügt.
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 = zeit; % 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;
% 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);
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 = current; % 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);
Diese Fehlermeldung besagt, dass du an irgend einer Stelle (welche du nicht nennst!) versuchst, unterschiedlich große Variablen miteinander zu verbinden.
ist current ein 1x2500 oder ein 2500x1 Array? Ich nehme mal an letzteres, welches auch dann zur Fehlermeldung passen würde, da der folgende Schritt nur mit einem Zeilen- nicht aber mit einem Spaltenarray funktioniert.
Wenn du 2500 Messwerte hast, solltest du in dem Code auch das N entsprechend anpassen. Die obere Codezeile mit dem Anfügen von Nullen würde ich weglassen.
Gruß DSP
Gast50
Gast
Beiträge: ---
Anmeldedatum: ---
Wohnort: ---
Version: ---
Verfasst am: 07.05.2016, 12:03
Titel:
Vielen Dank. Habe N auf 2500 angepasst und y=current' gemacht.
Die Diagramme werden jetzt alle geplottet.
Nun stellt sich bei mir die Frage, wie ich die Grundfrequenz und dessen Amplitude ablesen kann.
Ich weiß von vorneherein schon, dass mein Grundsignal z.B. eine Frequenz von 50Hz hat.
Soll ich dann im Plott "Amplitudengang nach Fensterung" die Amplitude bei 50Hz ablesen?
woher sollen wir wissen was du suchst bzw. was deine Grundschwingung ist? Wenn ein Peak bei 50 Hz ist, wäre das eben ein Signalanteil. Ob das nun die Grundfrequenz ist, kann man aus den gegebenen Infos nicht sagen. Am besten mal einen Plot des Amplitudenspektrums ins Forum stellen
Gruß DSP
Gast50
Gast
Beiträge: ---
Anmeldedatum: ---
Wohnort: ---
Version: ---
Verfasst am: 07.05.2016, 16:03
Titel:
Ich habe den aufgenommenen Strom und dessen Frequenzspektrum in der Datei angehängt.
Der Strom hat eine Frequenz von 30Hz.
Ich brauche jetzt die Grundschwingung um den Effektivwert den Grundschwingung bestimmen zu können.
Anschließend soll ich den Oberschwingungsgehalt sowie Grundschwingungsgehalt berechnen.
Bist du dir überhaupt sicher, dass die Abtastfrequenz deines Signals 8000 Hz ist? Oder hast du das einfach aus dem Skript übernommen? Wenn dein Signal mit einer anderen Abtastfrequenz aufgenommen wurde, ist das Ergebnis dann ohnehin nicht korrekt. Du musst schon die exakte Abtastfrequenz angeben.
Stimmt jetzt wo du es sagst.
Der Abtastfrequenz ist 2kHz.
Ich habe jetzt für die Abtastfrequenz 2000 eingegeben und dann wieder geplottet.
Das Ergebnis hänge ich mal an. Der Plot ist wieder mit f=30Hz.
Dieses mal ist bei dem Amplitudenspektrum nicht mehr viele Frequenzen zu sehen.
Ist meine Grundfrequenz nun die Frequenz die eine maximale Ampltuide hat?
Oder ist meine Grundfrequenz bei f=30Hz und von dort muss ich die Amplitude ablesen?
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