Ich will nun das Signal "b" in den Frequenzbereich transformieren und die Zeitverschiebung mit einem Totzeitglied korrigieren. Anschließend führe ich eine Rücktransformation durch um das Ergebnis zu kontrollieren:
Wenn ich nun ifftb und das Signal a plotte habe ich immernoch die zeitliche Verschiebung (0.01s). Nur die Amplitude des Signals "b" ist kleiner geworden. Außerdem ist im Signal ifftb immernoch ein komplexer Anteil. Dieser wird im plot vernachlässigt.
Habe ich etwas grundsätzlich falsch verstanden? Eine Verschiebung des Signals im Zeitbereich ist leider unerwünscht, dadurch würden die verwendbaren Signallängen verändert. Habe auch schon versucht mit der Funktion "tf" zu arbeiten. Aber da werde ich erst recht nicht schlau. Die ist nur ein vereinfachtes Beispiel. In der Praxis hätte ich z.B. zwei Sinus-Sweeps (Frequenz variiert), da hätte ich auch das zusätzliche Problem mit der zu verwendenden Frequenz, aber da könnt ich dann weitertesten, wenn dieses einfache Beispiel funktionieren würde... ich hoffe, dass ihr mir hierbei helfen könnt.
Ich grübel da auch schon seid gestern dran...aber selbst für dieses einfache Glied, reichen meine Regelungstechnik-Kenntnisse offenbar nicht aus.
Was aber definitiv geändert werden muss, dass man die Eulerform des Totzeitglieds benutzt. Addiere ich es jeweils zu Real- und Imaginärteil von Signal b, ändert sich nichts. Beim Multiplizieren, ändert sich sowohl Amplitude als auch Phase. Allerdings verstehe ich den Wert der Phasenverschiebung nicht. Das Signal b hat eine Verschiebung von 18 Grad, was einer Zeitdifferenz von 0.01s entspricht. Diese wäre ja auch meine Totzeit. Beim Totzeitglied ergibt sich eine Phasenverschiebung von
phi[rad] = - omega * Tt = 2pi * 5 Hz * Tt = -0.3142 in Grad sind das 18.0032.
Im Plot ist die Phasenverschiebung aber nun größer. Vielleicht kann ja mal ein Regelungstechniker Aufklärung leisten.
Code:
clearall;
Ts=0.0005;
t = 0:Ts:1;
% Signal
frequenz_a = 5;
frequenz_b = 5;
omega_a = 2*pi*frequenz_a;
omega_b = 2*pi*frequenz_b;
a = 4.5*sin(omega_a*t);
b = 4.5*sin(omega_b*t - 0.1*pi);
% Erstes Maximum von a und b ermitteln [Ca,Ia]=max(a);
[Cb,Ib]=max(b);
% Laufzeitverschiebung berechnen = Totzeit
delta_t = Ib*Ts - Ia*Ts
% Phasenverschiebung von b in Grad
phi = 360 * frequenz_b * delta_t
% Fouriertransformation
fftb = fft(b);
% Aufteilung in Real- und Imaginärteil
re = real(fftb);
im = imag(fftb);
% Totzeit: F(jw) = Kp * exp(-j w Tt) % Eulerscher Satz % F(jw) = Kp * ( cos(wTt) - jsin(wTt) ) % Totzeit zu Signal b
Kp = 1;%2.4322;
re = re .*(Kp * cos(omega_b*delta_t));
im = im .*(- Kp * sin(omega_b*delta_t));
fftc = complex(re,im);
ifftc = ifft(fftc);
hold on;plot(t,ifftc,'g*-');
ffta = fft(a);
% Amplitude von b
Mag_b = abs(fftb);
% Phase von a
phase_a = angle(ffta);
% Signal c aus a und b zusammensetzen
re = Mag_b .*cos(phase_a);
im = Mag_b .*sin(phase_a);
fftc = complex(re,im);
% Rücktransformation
ifftc = ifft(fftc);
hold on;plot(t,ifftc,'k--','LineWidth',2);
legend('Signal a', 'Signal b', 'Signal b + Totzeit', 'Signal c aus a und b')
Ich habe jetzt noch eine Lösung beigefügt, wie man die Phasenverschiebung anderweitig beseitigen kann. Allerdings macht das wenig Sinn, wenn die Signale a und b unterschiedliche Frequenzen haben.
Edit: Allerdings geht das natürlich auch wesentlich einfacher. Wenn man die Zeitdifferenz Tt von a und b kennt -> xa(t) = Kp * xe(t - Tt). Einfach den Zeitvektor ändern
Wenn ich mir allerdings das Bodedigramm des Totzeitglieds anschaue, macht das vorherige Ergebnis keinen Sinn. Die Verstärkung ist 1 und die Phase bei 5 Hz beträgt meine gewünschten 18 Grad.
Code:
f=0.1:0.1:10;
w = 2*pi.*f;
Kp = 1;
re = (Kp * cos(w*delta_t));
im = (- Kp * sin(w*delta_t));
Betrag = abs(complex(re,im));
phase = unwrap(angle(complex(re,im)));
figure(2) subplot(211) loglog(w,Betrag);grid on
subplot(212) semilogx(w,angledim(phase,'radians','degrees'));grid on
hold on;plot(2*pi*5,[0:-0.1:-20],'r','LineWidth',4)
Vermutlich ist also an der Rechnung noch etwas falsch. Oder man kann die Phasenverschiebung so nicht über das Totzeitglied realisieren.
Giuseppe
Gast
Beiträge: ---
Anmeldedatum: ---
Wohnort: ---
Version: ---
Verfasst am: 14.07.2011, 10:54
Titel:
So, hier meine Lösung, habe ich nun doch noch geschafft:
Totzeitglied wird angewandt, Zeitverschiebung ist bekannt. Wenn es nicht der Fall sein sollte könnte man diese durch eine Kreuzkorrelation detektieren.
Bemerkung: Lösung funktioniert aber löst nicht das eigentliche Problem:
Die Verschiebung klappt zwar, aber wenn man (wie in diesem Beispiel) ein Signal verarbeitet, welches nicht genau einem Vielfachen der Periode entspricht, dann sieht man einen Sprung im verschobenem Signal (in diesem Fall zum Schluss). Für ein aperiodisches Signal kann man also die Verschiebung auch im Zeitbereich durchführen... eigentlich auch logisch: das nach rechts verschobene Signal kann keine zusätzlichen Informationen (hier: am Ende des Signals) aus dem Hut zaubern: was an einer Seite im Zeitbereich "abgeschnitten" wird, wird auf der anderen Seite angefügt.
Danke für eure Antworten, da wurde ja richtig Zeit investiert, hoffe mich revanchieren zu können. Wenn Fragen aufkommen, werd ich die gerne beantworten
T = 0.0005;
Fs = 1/T;
L = length(b);
NFFT = fft(b)/L;
% Wenn man eine gerader Sample-Anzahl hat % entsteht im NFFT-Vektor bei L/2+1 ein reeller Peak, % der für die negative und positive Maximal-Frequenz % gilt. Dieser darf bei der Verschiebung nicht % verändert werden, er muss Reel bleiben!! % NFFT(1) darf nie verändert werden, unabhängig von % der Sampla-Anzahl. % Bei gerader Sample-Anzahl ist der Spiegelpunkt % bei (L+1)/2, hier muss der NFFT-Vektor ab NFFT(2) % komplett verschoben werden.
Wenn man Glieder der Regelungstechnik im Frequenzbereich anwenden will, kann man wohl immer so wie oben vorgehen. Glieder sind Frequenzabhängig und können daher nicht als konstante angesehen werden:
Jede Spektrallinie in NFFT ( =(fft(X) ) steht für eine Frequenz, diese muss zuerst festgestellt werden (Analoges Vorgehen wie bei der X-Achse für den Plot). Anschließend kann das Glied mit der jeweiligen Frequenz der Spektrallinie multipliziert werden.
Kann mir zwar nicht vorstellen, dass es nun eine neue Erkentnis hier im Forum ist, aber wenns öfter erwähnt wird, kann man es auch schneller finden
oho...ich Depp habe beim Totzeitglied immer nur das Omega des Signals eingesetzt, anstatt den Frequenzbereich der FFT zu nehmen. Da hätte ich aber spätestens beim Bodediagramm drauf kommen müssen
Auch das man Gleichanteil und Nyquistfrequenz nicht verändern darf, leuchtet ein. Danke für die Lösung! Das es bei nicht vollen Perioden im Messfenster nicht funktioniert, liegt ja aber nicht an der Rechnung/Idee mit dem Totzeitglied...sondern an dem ungenauen Spektrum "dank" Leakage der FFT.
Habe nun den Code auf das tatsächliche Problem angewandt:
Ich messe eine fremderregte Schwingung mit einem Beschleunigungssensor und mit einem optischen System. Das opt. System hinkt etwas hinterher (ca. 4ms).
Die Verschiebung hat nun wunderbar geklappt. Leider habe ich zum Schluss etwas ungewöhnliches (siehe Bild). Kann das vom Leakage kommen? Wenn ja, kann man das verhindern? Ich denke, dass es wohl nur mit einer höheren Samplingrate oder einer höheren Sample-Anzahl besser klappt, da dann das Frequenzband eine höhere Auflösung hat. Wär jetzt auch nicht schlimm. Muss ich das Ende halt trimmen... aber für mein Verständnis wäre das schon interessant. Man kann auch erkennen, dass das Original-Signal diesen Fehler nicht hat.
Probier doch mal das folgende Programm aus...ersetze das Testsignal durch dein verschobenes Messsignal des optischen Sensors. Das Programm stammt von hier aus einem anderen Thread, wo nach der Grundwelle eines Signals gesucht wurde.
Du kannst ja mal Schrittweise das k in Zeile 46 erhöhen, dann wird das Signal nicht nur aus der Grundwelle, sondern auch mit k-1 Oberwellen zusammengesetzt. Tritt der Fehler auch schon bei k = 1 auf?
signal = offset+(Ampl1*sin(2*pi*f1*zeit + phase1) + Ampl2*sin(2*pi*f2*zeit + phase2));
% reales Signal
%signal = y';
%zeit = x;
% 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));
phase = angle(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));
phase_(i) = phase(position_f(length(position_f)-i+1));
betrag_(i) = wert_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_;
phase = phase_;
betrag = betrag_;
%% Plot der Berechnungen figure(1) subplot(311) plot(zeit,signal,'r','LineWidth',2);
xlabel('Zeit [s]');
ylabel('Signal y(t)');
title('Zeitsignal','FontSize',18);
grid on
subplot(312) stem(f,abs(SIGNAL),'r','LineWidth',2,'Markersize',8);
xlabel('Frequenz [Hz]');
ylabel('Absolutwert');
title('Frequenzspektrum - Amplitude','FontSize',18);
grid on
subplot(313) plot(f,angledim(angle(SIGNAL),'radians','degrees'),'r.-','LineWidth',2,'Markersize',8);
xlabel('Frequenz [Hz]');
ylabel('Phase [Grad]');
title('Frequenzspektrum - Phase','FontSize',18);
grid on
figure(2) hold on
plot(zeit,signal,'r','LineWidth',2);
plot(zeit,signal_FFT,'b.-','LineWidth',1);
plot(zeit,signal_FFT_Phasenaenderung,'g--','LineWidth',4);
hold off
title('Ermittlung der Grundwelle des Signals','FontSize',18);
h=sprintf(['orig. Grundwelle mit Phase = ' num2str(angledim(phase(k),'radians','degrees'),'%3.1f') ' Grad'],'FontSize',14);
legend('original',h,'Grundwelle mit Phasenänderung');
disp('Zeitverschiebung');
delta_t = (phase(k)+(pi/2))/(2*pi*f_(k)) xlabel('Zeit [s]');
ylabel('Signal y(t)');
grid on
pause closeall
Danke für den Code. Leider bekomme ich bei der Ausführung den Fehler:
Error using ==> horzcat ; in Zeile 12
Kann mir den Fehler nicht wirklich erklären. Schaue ich mir morgen nochmal genauer an.
Wenn du diese Programm ohne Änderungen ausführst? Bei mir läuft es einwandfrei.
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