WICHTIG: Der Betrieb von goMatlab.de wird privat finanziert fortgesetzt. - Mehr Infos...

Mein MATLAB Forum - goMatlab.de

Mein MATLAB Forum

 
Gast > Registrieren       Autologin?   

Partner:




Forum
      Option
[Erweitert]
  • Diese Seite per Mail weiterempfehlen
     


Gehe zu:  
Neues Thema eröffnen Neue Antwort erstellen

Grundwelle eines Signals berechnen

 

DSP
Forum-Meister

Forum-Meister



Beiträge: 2.117
Anmeldedatum: 28.02.11
Wohnort: ---
Version: R2014b
     Beitrag Verfasst am: 05.07.2011, 11:50     Titel:
  Antworten mit Zitat      
Das ist ja alles richtig und mir auch bekannt. Meine Info über die Phase bezog sich auch nicht auf die Berechnung deines Zeitsignals sondern lediglich darauf, dass er nicht diesen Phasenwert von

Code:
Phase = angle(SIGNAL) % in radians


zum Schieben verwenden kann Wink. Er bekommt zwar eine Phasenverschiebung, weiß dann aber nicht ob die Grundschwingung ein Sinus oder Cosinus war (mal vom Bild abgesehen).

Ich habe es aber gerade nochmal getestet und so geht es doch...egal ob Sinus oder Cosinus als Grundschwingung. Beim Cosinus gibt es dann aber eine Phasenverschiebung von 90 Grad. Aber die Grundwelle beginnt immer beim Nulldurchlauf.


Code:

[phase] = sort(angle(SIGNAL));
.
.
.
for i=1:1:length(zeit)
     summand     = 0;
     for m=1:1:k
         summand     = a(m)*cos(2*pi*f_(m)*zeit(i) + phase(m)) + b(m)*sin(2*pi*f_(m)*zeit(i) + phase(m)) + summand;
     end;
     signal_FFT(i)       = a_0 + summand;
end;
 
Private Nachricht senden Benutzer-Profile anzeigen


Idefix_1024
Forum-Century

Forum-Century


Beiträge: 230
Anmeldedatum: 16.10.08
Wohnort: ---
Version: ---
     Beitrag Verfasst am: 05.07.2011, 12:03     Titel:
  Antworten mit Zitat      
ich verstehe weiterhin nicht was das Ziel mit dieser Phase sein soll...

Wenn ich den Code um die Phase wie von Dir vorgeschlagen erweitere ergeben sich verschobene Zeitverläufe, aber das Signal beginnt auch nicht bei null. Vielleicht kannst Du mal das gesamte Beispiel posten?
Private Nachricht senden Benutzer-Profile anzeigen
 
DSP
Forum-Meister

Forum-Meister



Beiträge: 2.117
Anmeldedatum: 28.02.11
Wohnort: ---
Version: R2014b
     Beitrag Verfasst am: 05.07.2011, 12:44     Titel:
  Antworten mit Zitat      
Die Phase gibt doch genau an, wie das Signal verschoben ist. Wenn ich sie kenne, kann ich das Signal auch so manipulieren, dass die Phasenverschiebung weg ist. Ob es Sinn macht, ist ja eine ganze andere Frage.

Es geht so aber doch nicht ganz. Ich habe es nur kurz getestet und mich auf Grund des Offsets des Testsignals täuschen lassen. Wenn man den Offset wie hier weglässt, treffe ich nicht exakt den Nulldurchgang. Obwohl meine Phasenverschiebung der Grundschwingung f1 60 Grad ist, beträgt sie mit dem Befehl angle() über 60 Grad Question Außerdem muss die Berechnung doch für Sinus und Cosinus als Testsignal unterschiedlich sein(wenn man die Phase manipulieren will). Diese Info hat man ja aber in der Regel nicht.

Code:

%% Signal
 dt          = 1/8000;
 zeit        = 0:dt:0.1;
 f1          = 100;
 f2          = 500;
 Ampl1       = 10;
 Ampl2       = 5;
 offset      = 0;
 phase1      = (pi/3);
 phase2      = 0;

 signal      = offset+(Ampl1*cos(2*pi*f1*zeit + phase1) + Ampl2*cos(2*pi*f2*zeit + phase2));


% 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 position_f2]   = sort(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_f2(length(position_f2)-i+1));
end

% Signal auf Grund der FFT-Werte erzeugen bis zur k-ten Harmonischen
 k           = 1;
 a           = a_;
 b           = b_;
 phase       = phase_;

% Offset separat behandeln!
 gleichanteil = 0;
if find(f_==0)
     a_0             = a(find(f_==0));
     a(find(f_==0))  = [];
     gleichanteil    = 1;
end;

for i=1:1:length(zeit)
     summand     = 0;
     for m=1:1:k
         summand     = a(m)*cos(2*pi*f_(m)*zeit(i) + phase(m)) + b(m)*sin(2*pi*f_(m)*zeit(i) + phase(m)) + summand;
     end;
     signal_FFT(i)       = a_0 + summand;
end;


% Signal-Berechnung auf Grund eines Filters
 Anzahl_mw = 64;
 b(1:Anzahl_mw) = 1/64;
 signal_filt = filter(b,1,signal);
 

 %% Plot der Berechnungen
figure(1)
subplot(311)
plot(zeit,signal,'r','LineWidth',2);
xlabel('Zeit [s]');
ylabel('Signal');
grid on
subplot(312)
plot(f,abs(SIGNAL),'rx','LineWidth',2,'Markersize',8);
xlabel('Frequenz [Hz]');
ylabel('Absolutwert');
grid on
subplot(313)
plot(f,angledim(angle(SIGNAL),'radians','degrees'),'r.-','LineWidth',2,'Markersize',8);
xlabel('Frequenz [Hz]');
ylabel('Phase [Grad]');
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
 
Private Nachricht senden Benutzer-Profile anzeigen
 
Idefix_1024
Forum-Century

Forum-Century


Beiträge: 230
Anmeldedatum: 16.10.08
Wohnort: ---
Version: ---
     Beitrag Verfasst am: 05.07.2011, 13:49     Titel:
  Antworten mit Zitat      
wenn Du den Zeitvektor auf 0,01 begrenzt erhältst Du "genau" eine Periode der Schwingung. Dann siehst Du in Deinem Phasen-Plot bei 100 Hz eine Phase von 1,07 (= pi/3 mit Unschärfe) und bei 500 Hz hast Du eine Phase von Null (wieder nicht exakt aber ziemlich gut).

Wenn Du aber nun das zeitliche Signal über einen längeren Zeitraum betrachtest kannst Du das so nicht mehr herauslesen. Dann verschmiert sich das Ergebnis genau so wie es bei Deinem Plot zu sehen ist.

aso, zur Info, Dein Plot ist nicht in Grad sondern in Radians

ist also alles wie es sein soll...

Grüße,
Idefix_1024
Private Nachricht senden Benutzer-Profile anzeigen
 
DSP
Forum-Meister

Forum-Meister



Beiträge: 2.117
Anmeldedatum: 28.02.11
Wohnort: ---
Version: R2014b
     Beitrag Verfasst am: 05.07.2011, 14:16     Titel:
  Antworten mit Zitat      
Das mit der Unschärfe wusste ich nicht...wieder was dazu gelernt. Wenn das Testsignal jetzt aber ein Sinus ist, geht es mit dem Addieren der Phase so nicht mehr. Dann habe ich keinen Nulldurchgang bei t=0.

Ich habe hier aber wohl ohnehin einen Denkfehler drin, da in den Koeff. a und b ja schon die Info über die Phase steckt. Also müsste man eigentlich a und b manipulieren. Aber es gibt ja auch folgende Beziehung...

Code:
a*cos(x)+b*sin(x) = betrag*cos(x+phase)


Wenn mich nun nur die Amplitude und Frequenz der Grundschwinung interessiert, kann ich ja auch die Phase weglassen. Da er ja aber bei t=0 einen Nulldurchgang haben möchte, müsste ich die Phase = pi/2 setzen...da es ja ein cos ist.

Warum sollten das nicht Grad sein? Ich rechne doch mit angledim in Grad um, da angle mir den Winkel im Bogenmaß liefert.
Private Nachricht senden Benutzer-Profile anzeigen
 
DSP
Forum-Meister

Forum-Meister



Beiträge: 2.117
Anmeldedatum: 28.02.11
Wohnort: ---
Version: R2014b
     Beitrag Verfasst am: 05.07.2011, 14:44     Titel:
  Antworten mit Zitat      
Nun geht es bei cos und sin als Testsignal (auch mit längerer Messzeit). Ich hatte beim sortieren von Betrag und Phase noch einen Fehler. Jetzt kommt allerdings immer ein Sinus raus...

Code:

%% Signal
clear all;
 dt          = 1/8000;
 zeit        = 0:dt:0.1;
 f1          = 100;
 f2          = 500;
 Ampl1       = 10;
 Ampl2       = 5;
 offset      = 0;
 phase1      = (pi/3);
 phase2      = 0;

 signal      = offset+(Ampl1*cos(2*pi*f1*zeit + phase1) + Ampl2*cos(2*pi*f2*zeit + phase2));


% 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 position_f2]   = sort(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(length(position_f2)-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_;

% Offset separat behandeln!
 gleichanteil = 0;
if find(f_==0)
     a_0             = a(find(f_==0));
     a(find(f_==0))  = [];
     gleichanteil    = 1;
end;

%phase(1)=1.04719755119660;
for i=1:1:length(zeit)
     summand     = 0;
     for m=1:1:k
         summand     = betrag(m)*cos(2*pi*f_(m)*zeit(i) - (pi/2)) + summand;
     end;
     signal_FFT(i)       = a_0 + summand;
end;


% Signal-Berechnung auf Grund eines Filters
 Anzahl_mw = 64;
 b(1:Anzahl_mw) = 1/64;
 %h_sinc = window_sinc_filter(100,400,8000,1) ;
 signal_filt = filter(b,1,signal);
 

 %% Plot der Berechnungen
figure(1)
subplot(311)
plot(zeit,signal,'r','LineWidth',2);
xlabel('Zeit [s]');
ylabel('Signal');
grid on
subplot(312)
plot(f,abs(SIGNAL),'rx','LineWidth',2,'Markersize',8);
xlabel('Frequenz [Hz]');
ylabel('Absolutwert');
grid on
subplot(313)
plot(f,angledim(angle(SIGNAL),'radians','degrees'),'r.-','LineWidth',2,'Markersize',8);
xlabel('Frequenz [Hz]');
ylabel('Phase [Grad]');
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);
%plot(zeit(1:length(signal)-32),signal_filt(1,33:length(signal)),'k--','LineWidth',2);
hold off
legend('original','FFT','Filter',0)
xlabel('Zeit [s]');
ylabel('Signal');
grid on
pause
close all
 
Private Nachricht senden Benutzer-Profile anzeigen
 
Rickson
Themenstarter

Forum-Newbie

Forum-Newbie


Beiträge: 8
Anmeldedatum: 29.06.11
Wohnort: ---
Version: ---
     Beitrag Verfasst am: 08.07.2011, 23:50     Titel:
  Antworten mit Zitat      
Hallo DSP, hallo Idefix_1024

Der von DSP gepostete Code scheint sehr gut zu funktionieren Very Happy

Mittlerweile beschäftigt mich noch eine andere Fragestellung.
Hierfür habe ich euch mal das .mat-File des verrauschten, sinusförmigen Messsignals angehängt.

Es ist im Plot deutlich zu erkennen, dass das genannte Signal eine Phasenverschiebung besitzt.

Ich frage mich nun, ob diese Phasenverschiebung durch horizontales Verschieben des Signals auf der Zeitachse eliminiert werden kann.

Hier stößt man jedoch auf die folgenden Probleme:
1.) Wie ermittelt man eigentlich die vorhandene Phasenverschiebung des Signals, die es zu eliminieren gilt?
2.) Wie kann ich das Signal anschließend so auf der Zeitachse verschieben, dass ein Signal entsteht, welches:
a) bei x=0 seinen positiven Nulldurchgang hat und ansonsten unverändert bleibt?
b) dem Signal auch a) um 90° voreilt (Cosinus)?

Ich vermute, dass sich diese Aufgabe ebenfalls im Frequenzbereich lösen lässt indem man z.B. nach der fft jede Komponente mit exp(i*p) multipliziert (wobei p die in 1.) ermittelte Phasenverschiebung ist).

Die Implementierung will mir aber bisher nicht gelingen.
Habt ihr vielleicht eine gute Idee?

Ein schönen Abend!

matlab.zip
 Beschreibung:

Download
 Dateiname:  matlab.zip
 Dateigröße:  4.18 KB
 Heruntergeladen:  625 mal
Private Nachricht senden Benutzer-Profile anzeigen
 
Rickson
Themenstarter

Forum-Newbie

Forum-Newbie


Beiträge: 8
Anmeldedatum: 29.06.11
Wohnort: ---
Version: ---
     Beitrag Verfasst am: 09.07.2011, 09:44     Titel:
  Antworten mit Zitat      
Da hat sich noch ein Rechtschreibfehler eingeschlichen. Es sollte heißen:

"b) dem Signal AUS a) um 90° voreilt (Cosinus)? "

Embarassed
Private Nachricht senden Benutzer-Profile anzeigen
 
DSP
Forum-Meister

Forum-Meister



Beiträge: 2.117
Anmeldedatum: 28.02.11
Wohnort: ---
Version: R2014b
     Beitrag Verfasst am: 09.07.2011, 17:00     Titel:
  Antworten mit Zitat      
Ich hatte doch noch einen Fehler bei der Phase, der aber jetzt behoben ist. Ehrlich gesagt verstehe ich deine Fragen nicht. Ich habe nun dein beigefügtes Signal getestet und die Plots entsprechend angepasst. Die Phasenverschiebung der Grundwelle deines Zeisignals wird in der Legende angezeigt. Bei diesem Signal betrifft die Phasenverschiebung -140 Grad.

Da ich für die Änderungen aber die Phase weglasse erhältst du doch, ein Signal mit positivem Nulldurchgang bei t=0...was nur ein Sinus sein kann.
Warum willst du jetzt noch schieben? Es ist doch auch gegenüber dem Original in Amplitude und Frequenz identisch...nur die Phase ändere ich.

Die Phasenverschiebung berechnet man aus den Fourierkoeffizienten a und b = Real- und Imaginärteil so.

phi = arctan(b/a) -> phi in [rad]

in Matlab nutzt man einfach die Funtkion angle() oder um es wie hier angegeben zu berechnen, atan2(). Ich gebe dir im Command Window auch noch die Zeitverschiebung delta_t = (phi + 90 Grad)/freq. aus. Das Signal muss um diesen Wert nach links auf der Zeitachse verschoben werden, um keine Phasenverschiebung zu haben und somit deine Forderung erfüllt zu haben.

Teste doch einfach nochmal mit dem Testsignal rum, damit es klar wird. Der Cosinus eilt dem Sinus um 90 Grad voraus. Nehmen wir mein Psp. mit den pi/3 = 60 Grad Phasenverschiebung bei f1 an. Ist das Testsignal ein Sinus, erhältst du theoretisch -30 Grad Phase (real auf Grund von Ungenauigkeiten nicht ganz) aus der oben benannten Berechnung. Wenn du zu den -30 nun 90 Grad Verschiebung zum Cosinuns dazurechnest, erhältst du wieder die ursprunglichen 60 Grad deines Testsignals. Ist das Testsignal degegen ein Cosinus, wirst du gleich aus der Berechnung die 60 Grad erhalten.

Code:

% Testsignal
 %dt          = 1/10000;
 %zeit        = 0:dt:0.01;
 %f1          = 100;
 %f2          = 500;
 %Ampl1       = 10;
 %Ampl2       = 5;
 %offset      = 0;
 %phase1      = (pi/3);
 %phase2      = 0;

 %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_;

% Offset separat behandeln!
 gleichanteil = 0;
if find(f_==0)
     a_0             = a(find(f_==0));
     a(find(f_==0))  = [];
     gleichanteil    = 1;
end;

for i=1:1:length(zeit)
     summand     = 0;
     summand2    = 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;
         summand     = betrag(m)*cos(2*pi*f_(m)*zeit(i) + phase(m)) + summand;
         summand2    = betrag(m)*cos(2*pi*f_(m)*zeit(i) - (pi/2)) + summand2;
     end;
     signal_FFT(i)                 = a_0 + summand;
     signal_FFT_Phasenaenderung(i) = a_0 + summand2;
end;

%% 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)
plot(f,abs(SIGNAL),'rx','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');
xlabel('Zeit [s]');
ylabel('Signal y(t)');
grid on

disp('Zeitverschiebung');
delta_t = (phase(k)+(pi/2))/(2*pi*f_(k))

pause
close all
 
Private Nachricht senden Benutzer-Profile anzeigen
 
Neues Thema eröffnen Neue Antwort erstellen

Gehe zu Seite Zurück  1, 2, 3

Einstellungen und Berechtigungen
Beiträge der letzten Zeit anzeigen:

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
.





 Impressum  | Nutzungsbedingungen  | Datenschutz | FAQ | goMatlab RSS Button RSS

Hosted by:


Copyright © 2007 - 2024 goMatlab.de | Dies ist keine offizielle Website der Firma The Mathworks

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.