äh, da hab ich das wohl gestern abend falsch gelesen, mit dem Signal wiederherstellen. War der festen Überzeugung, den TE interessiert nur die Frequenz...
okay, also, f hat man, die Nullstellen als Median aller Nulldurchgänge auch (damit die Phasenlage). Damit sollte es möglich sein, einen sinus-fit selber zu schreiben, oder?
Erste Näherung durch die absoluten Maxima gibts auch,
von da aus die Amplitude verringern,
bis die Abweichung Signal -- Fit wieder größer wird,
fertig.
[edit:] noch ne idee (allgemein) wäre einen gleitenden Mittelwert von rechts und links drüber zu legen, und aus dem arithmetischen Mittel der beiden Gleitenden dann eine gute Näherung zu haben.
_________________
Ich hasse es wenn die Leute Fragen stellen, man dann versucht sich Mühe zu geben, und diejenigen ihren Thread nie wieder besuchen...
na da stimmt mir wohl jemand zu nur weil ich "viel" code gepostet habe...
Mein code war so ja gar nicht lauffähig! Ist mir gerade aufgefallen! Tut mir leid, mir war nicht aufgefallen, dass ich eine Funktion aufgerufen habe, die Ihr gar nicht haben könnt...
Deshalb nun der Quellcode, der für sich lauffähig sein sollte!
Code:
%% Signal
dt = 1/8000;
zeit = 0:dt:0.1;
f1 = 100;
f2 = 500;
Ampl1 = 10;
Ampl2 = 5;
offset = 2;
signal = offset+(Ampl1*cos(2*pi*f1*zeit)+Ampl2*cos(2*pi*f2*zeit));
% FFT Berechnen
Ts = diff(zeit(1:2)); % Abtastzeit
N = length(zeit); % Länge des Datenvektors
f = [0:floor((N-1)/2)] / (N*Ts); % Frequenzvektor für Spektrum-Plot % Fouriertransformierte
SIGNAL = fft(signal); % Fouriertransformation
SIGNAL = SIGNAL / N; % Normierung
SIGNAL = [SIGNAL(1)2*SIGNAL(2:floor((N-1) / 2) + 1)]; % begrenzen auf f_max
% Fourierkoeffizienten a und b bestimmen for i=1:1:length(SIGNAL)
a(i) = real(SIGNAL(i));
b(i) = -imag(SIGNAL(i));
end;
% Frequenzvektor sortieren so dass später das Signal aus den Hauptfrequenzanteilen % wieder erzeugt bzw. approximiert werden kann [wert_f position_f] = sort(abs(SIGNAL));
for i=1:1:length(f)
f_(i) = f(position_f(length(position_f)-i+1));
a_(i) = a(position_f(length(position_f)-i+1));
b_(i) = b(position_f(length(position_f)-i+1));
end
% Signal auf Grund der FFT-Werte erzeugen bis zur k-ten Harmonischen
k = 1;
a = a_;
b = b_;
for i=1:1:length(zeit)
summand = 0;
if gleichanteil == 1 for m=1:1:k
summand = a(m)*cos(2*pi*f_(m)*zeit(i)) + b(m)*sin(2*pi*f_(m)*zeit(i)) + summand;
end;
end;
if gleichanteil == 0 for m=1:1:k
summand = a(m)*cos(2*pi*f_(m)*zeit(i)) + b(m)*sin(2*pi*f_(m)*zeit(i)) + summand;
end;
end;
signal_FFT(i) = a_0 + summand;
end;
% Signal-Berechnung auf Grund eines Filters
Anzahl_mw = 32;
b(1:Anzahl_mw) = 1/32;
signal_filt = filter(b,1,signal);
%% Plot der Berechnungen figure(1) subplot(211) plot(zeit,signal,'r','LineWidth',2);
xlabel('Zeit [s]');
ylabel('Signal');
grid on
subplot(212) plot(f,abs(SIGNAL),'rx','LineWidth',2,'Markersize',8);
xlabel('Frequenz [Hz]');
ylabel('Absolutwert');
grid on
figure(2) hold on
plot(zeit,signal,'r','LineWidth',2);
plot(zeit,signal_FFT,'b','LineWidth',2);
plot(zeit,signal_filt,'m','LineWidth',2);
hold off
legend('original','FFT','Filter',0) xlabel('Zeit [s]');
ylabel('Signal');
grid on
und wenn man den Code nun ausführt, dann merkt man, dass das mit dem Filter tatsächlich seine Tücken hat...
dass es aber einfacher zu implementieren ist, ist keine Frage ;-)
Ich habe den Code irgendwann geschrieben, weil mich interessiert hat, wie ein Signal aussieht wenn es nur die Grundwelle und die x-te Harmonische beinhaltet. Vielleicht nutzt es hier ja auch sonst noch jemandem!
Ich brauche dir sicherlich nicht nach dem Mund reden...das dein Code nicht funktioniert wegen dem fehlenden Code für die Funktion, habe ich gesehen. Ich hatte das Filter bei mit selber schon getestet und meinen Fehler sehen können. Das Filter hatte ich einfach als Bsp. mit angegeben, ohne es vorher überhaupt getestet zu haben oder an die kurze Messdauer zu denken. Als ich dann meinen Fehler gepostet habe, war dein neuer Post/ Code noch gar nicht zu sehen
Edit: das Messfenster ist aber immer noch größer als die Grundperiode von f1
Ein längerer Zeitraum (bzw. mehrere Perioden) führt bei einer FFT allerdings eigentlich zu genaueren Ergebnissen.
Hier der Plot in dem man den Fehler aller Verfahren sehen kann. Der Verlauf "ideal" würde sich ergeben, wenn man das Signal wirklich perfekt treffen würde.
Aso, nochmals Entschuldigung an DSP! Das Filter wäre natürlich deutlich besser als hier gezeigt, wenn man das Signal sinusförmig und nicht cosinusförmig wählen würde! Dann hätte man ein kleineres Einschwingproblem (siehe zweites Bild).
Ich bin sicherlich kein Lehrmeister, aber ein Anfänger auch nicht. Das Filter passt hier halt wirklich nicht, wenn nur eine Periode vorhanden ist. Wären es mehr, konnte man ja auch die Phasenverschiebung des Filters einfach entfernen. Ich hatte das Filter lediglich als möglichen Ansatz gepostet, zu mal das MW hier auch mit mehreren Perioden "überfordert" ist. Hätte er deinen schönen Code nicht, wäre es für ihn ohne Kenntnisse über die FFT sicherlich schwer geworden so etwas zu erstellen.
Wenn das Signal aber keine kompl. Periode(n) enthält, geht deine Funktion wegen Leakage aber auch nicht perfekt. Bei einer Periode kommt man dann auch mit Fensterung nicht wirklich weiter. Hier ist dann das Filter wieder im Vorteil.
Hallo!
Der Beitrag ist ja auf reges Interesse gestossen
Hätte nun noch eine weitere Frage zum Code von Idefix_1024:
Code:
%% Signal
dt = 1/8000;
zeit = 0:dt:0.1;
f1 = 100;
f2 = 500;
Ampl1 = 10;
Ampl2 = 5;
offset = 2;
signal = offset+(Ampl1*cos(2*pi*f1*zeit)+Ampl2*cos(2*pi*f2*zeit));
Ts = diff(zeit(1:2)); % Abtastzeit
N = length(zeit); % Länge des Datenvektors
f = [0:floor((N-1)/2)] / (N*Ts); % Frequenzvektor für Spektrum-Plot % Fouriertransformierte
SIGNAL = fft(signal); % Fouriertransformation
SIGNAL = SIGNAL / N; % Normierung
SIGNAL = [SIGNAL(1)2*SIGNAL(2:floor((N-1) / 2) + 1)]; % begrenzen auf f_max
ergebnis = [SIGNAL; f];
%% Plot der Berechnungen figure() subplot(211) plot(zeit,signal,'r','LineWidth',2);
xlabel('Zeit [s]');
ylabel('Signal');
grid on
subplot(212) plot(f,abs(SIGNAL),'rx','LineWidth',2,'Markersize',8);
xlabel('Frequenz [Hz]');
ylabel('Absolutwert');
grid on
Kann man (durch weiteren Code) aus dem berechneten Amplitudengang (also ohne auf den Plot zu schauen) alle erforderliche Kenngrössen extrahieren um die Grundwelle des Signals plotten zu können?
Ich habe die Koeffizienten a_ und b_ sortiert, so dass man bequem beliebige Anteile im Spektrum weglassen und berücksichtigen kann.
Dazu muss man den Faktor k abändern (derzeit auf 1 für "nur Grundwelle berücksichtigen").
An dieser Stelle muss ich aber sagen, dass man das Resultat besser nochmal im plot kritisch betrachten sollte, da der Code nicht wirklich getestet wurde ;-)
Gerade die Berücksichtigung eines eventuellen Offsets ist recht spontan eingebaut worden!
Hallo!
Vielen Dank nochmals für eure Beiräge. Die Grundwelle konnte ich mittlerweile erfolgreich berechnen. Auch habe ich mich nun etwas
tiefer in den FFT-Algorithmus eingelesen.
Leider beginnt die berechnete Grundwelle jedoch bei t=0 bei einem Wert ungleich Null was vermutlich damit zu tun hat, das das
verrauschte Signal (welches ja die Grundlage für die FFT bildete) ebenfalls nicht exakt bei Null beginnt.
Es besitzt vielmehr eine Phasenverschiebung und liegt in folgender Form vor:
[0.02; 0.07; 0.12; 0.19; 0.23 ; ... ; -0.08; -0.065; -0.02; -0.015].
Gibt es irgendeine Möglichkeit, das Signal so zu modifizieren. dass es bei y(t=0) = 0 beginnt und bei y(t=T) = 0 endet?
Mit der Phasenverschiebung kannst Du nun das Signal so "hinschieben" dass es Deinen Wünschen entspricht. Wenn man weit genug schiebt sieht es dann auch irgendwann mal aus wie ein cosinus Signal.
Man sollte sich aber schon im Klaren sein, dass man damit das Signal nachträglich, also künstlich verfälscht!!
Ich habe noch eine Frage an Idefix...kann es sein, dass du nach copy und paste was vergessen hast? Für mich sind die Ausdrücke in den if Abfragen identisch
Code:
for i=1:1:length(zeit)
summand = 0;
if gleichanteil == 1 for m=1:1:k
summand = a(m)*cos(2*pi*f_(m)*zeit(i)) + b(m)*sin(2*pi*f_(m)*zeit(i)) + summand;
end;
end;
if gleichanteil == 0 for m=1:1:k
summand = a(m)*cos(2*pi*f_(m)*zeit(i)) + b(m)*sin(2*pi*f_(m)*zeit(i)) + summand;
end;
end;
signal_FFT(i) = a_0 + summand;
end;
Allerdings ist das mit der berechneten Phase auch nicht ganz so einfach...du kannst sie nicht einfach, wie Idefix hier beschrieben hat, dazu addieren. Du willst ja immer einen Sinuns ohne Phasenverschiebung, so dass y(t=0)=0 ist.
Nehmen wir nun mal als Bsp. einen Cosinus Signal mit pi/3 = 60 Grad Phasenverschiebung als deine Grundschwingung in deinem gemessenen Signal an. Das es ein Cosinuns ist, kannst du aber nicht wissen. Die Phasenverschiebung von 60 Grad wirst du aus dem Spektrum der FFT erhalten. Der Cosinus ist aber 90 Grad Phasenverschoben zum Sinus. Diese Info bekommst du aus der Phase nicht und du weißt ebenfalls nicht, ob du einen Sinus oder Cosinus als Grundwelle hast.
Deine Ausführungen zur Phase verstehe ich ehrlich gesagt nicht wirklich.
Die FFT liefert zwei Koeffizienten, die bei mir a und b heißen. Der Koeffizient a beschreibt den cosinus Anteil und der Koeffizient b beschreibt den sinus Anteil. Genauer gesagt ergibt sich der cosinus Anteil aus dem Realteil und der sinus Anteil aus dem Imaginärteil des FFT Ergebnisses.
Das klappt IMMER ganz gleich ob das Ausgangssignal ein cos oder sin ist. Du kannst in meinem Beispiel-Code einfach aus sinus einen cos machen in den ersten Zeilen... also
Code:
signal = offset+(Ampl1*cos(2*pi*f1*zeit)+Ampl2*cos(2*pi*f2*zeit));
signal_ideal = offset+Ampl1*cos(2*pi*f1*zeit);
Jedes Signal kann laut Fourier durch eine Überlagerung aus Sinus und Kosinus dargestellt werden. Das klappt auch bei zB Rechtecken...
also wie gesagt... evtl. verstehe ich Dich einfach falsch.
Natürlich darf man ein Signal nicht einfach mal eben verschieben, weil es dann nicht mehr das gemessene Signal ist. Um aber genau eine Periode darzustellen kann man durchaus das Signal passend schieben oder eben ausschneiden. Warum muss das Signal denn genau bei t=0 auch null sein in diesem Fall? Vielleicht hilft das bei der Aufklärung...
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