Phasenverschiebung mit Cross Correlation in Matlab

[Owen S. (senmeis)]

Servus,

mit dem folgenden Code wird die Cross Correlation von zwei 
Sinus-Signalen berechnet, die eine Phasenverschiebung von pi/4 haben. 
Theoretisch kann diese pi/4 Phasenverschiebung anhand der Scheitelwerte 
der Cross Correlation bestimmt werden. Wie macht man weiter mit diesen 
Cross Correlation Werten in 'y' in diesem Beispiel?

hcorr = dsp.Crosscorrelator;
 t=[0:0.001:1];
 x1=sin(2*pi*2*t);
 x2=sin(2*pi*2*t – pi/4);
 y=step(hcorr,x2,x1); %computes cross-correlation of x1 and x2
 figure,plot(t,x1,'b',t, x2, 'g'); 
 legend('Input signal 1',' Input signal 2') 
 figure, plot(y); title('Correlated output')

Gruss
Senmeis

[Detlef _. (detlef_a)]

"Wie macht man weiter mit diesen Cross Correlation Werten in 'y' in 
diesem Beispiel?"

Was soll denn diese Frage bedeuten? Wenn Du mal zielgerichtete Fragen 
stellen würdest könnte man mal was antworten.

http://www.catb.org/esr/faqs/smart-questions.html

Cheers
Detlef

[Owen S. (senmeis)]

Sorry, ich habe einfach zu grob geschrieben.

Ich wollte eigentlich wissen, wie die Phasenverschiebung aus der 
Kreuzkorrelation berechnet werden soll. Diese Verschiebung ist irgendwie 
mit Stellen der Scheitelwerte verbunden.

Senmeis

[Frank (Gast)]

1)Kreuzkorrelation deiner Signale
2)Detektion um wieviele Abtastschritte das Maximum der KK verschoben 
ist.
3)aus der Abtastrate und Anzahl der Verschiebung deltaT berechnen.
4)auf den Phasenwinkel zurück rechnen.

[Owen S. (senmeis)]

OK, in diesem Beispiel:

t_sample = 0.001;
period = 0.5;
local_maximum = 157; % Index des ersten Maximums
deltaT = t_sample*local_maximum

Die Phase wird wie folgt berechnet:

deltaPhi = deltaT*2*pi/period;  % Ergebnis: 1.9729

Das stimmt aber nicht mit pi/4 = 0.7854.

Wo habe ich Fehler gemacht?

Senmeis

[Frank (Gast)]

Habs mal kurz probiert, sollte einigermaßen passen.

clear all;
close all;

%=================
%   Einstellungen
%=================
n=10;               %Anzahl Perioden
f=50;               %Signalfrequenz
fs=f*1000;         %Abtastfrequenz
phiGRAD = 50;       %Gewuenschte Phasenverschiebung in GRAD

%=================
%   Konstanten
%=================
radToGrad = 180/pi;
gradToRad = 1/radToGrad;
%=================
%   Berechnungen
%=================
T=1/f;              %Periodendauer
ts=1/fs;            %Abtastdauer
t=0:ts:(T*n - ts);
phiRAD = phiGRAD * gradToRad;

x1=sin(2*pi*f*t);
x2=sin(2*pi*f*t - phiRAD);

corr = conv(x1,x2);             %Faltung
tCorr = [fliplr(-t(2:end)) t];  %Zeitvektor fuer Faltung

h=figure('Name','Kreuzkorreltation zweier Signale');
maxCorrAmp = max(abs(corr));                %Absoluter Maximalwert der Kreuzkorrelation
[i,j,v] = find(abs(corr) == maxCorrAmp);    %Verschiebung suchen
deltaT = tCorr(j+1);                          %Verschiebung berechnen

phiCalcRad = 2*pi*f*deltaT;                  %aus KK berechnete Phase
phiCalcGrad = phiCalcRad * radToGrad;        %in Grad
phiErr = abs(phiCalcGrad-phiGRAD);            %Phasen-Fehler

%=================
%   Plots
%=================
subplot(2,1,1);
plot(t,x1,t, x2); 
legend('Input signal 1',' Input signal 2');
grid on;

subplot(2,1,2);
plot(tCorr,corr);
text(tCorr(j), corr(j), strcat('\leftarrow \Deltat = ', num2str(deltaT)), 'FontSize', 10);
legend('Crosscorrelation');
grid on;

%=================
%   Cmd
%=================
fprintf('delta_T =  %f s\n',deltaT);
fprintf('phiCalcGrad =  %f °\n',phiCalcGrad);
fprintf('phiErr =  %f °\n',phiErr);

[Owen S. (senmeis)]

Danke. Dein Code funktioniert.

Warum sollte conv() statt dsp.Crosscorrelator verwendet werden? In 
Deinem Code wird Faltung allerdings als KK kenngezeichnet.

Senmeis

[Frank (Gast)]

Owen Senmeis schrieb:
> Warum sollte conv() statt dsp.Crosscorrelator verwendet werden

dsp.Crosscorrelator ist in der DSP System Toolbox.
xcorr() in der Signal Processing Toolbox.
Habe ich beides nicht.
Mehr zum Thema:
http://dsp.stackexchange.com/questions/2654/what-is-the-difference-between-convolution-and-cross-correlation

Es müsste also auch xcorr() gehen, evtl muss man dann aber Anpassungen 
machen. Es sind sowieso wahrscheinlich nicht alle Fälle abgefangen!

Owen Senmeis schrieb:
> Deinem Code wird Faltung allerdings als KK kenngezeichnet.

Ja, das ist Formal nicht richtig...

[Owen S. (senmeis)]

Ich habe Abtastrate verringert und Fehler ändern sich dramatisch:

fs=f*1000 => phiErr =0.04
fs=f*100 => phiErr =7.6
fs=f*10 => phiErr =58

Gibt’s irgendwelche Optimierungsmöglichkeiten, Fehler für niedrige 
Abtastraten zu minimieren?

Senmeis

[Frank (Gast)]

Owen Senmeis schrieb:
> Ich habe Abtastrate verringert und Fehler ändern sich dramatisch:

Klar, Du hast ja auch weniger Informationen.

Owen Senmeis schrieb:
> Gibt’s irgendwelche Optimierungsmöglichkeiten, Fehler für niedrige
> Abtastraten zu minimieren?

Mein erster Ansatz wäre eine Interpolation der Zeitsignale bevor man die 
KK rechnet, z.B. mit einem Spline.

Evtl. könnte man auch mit einer Variation der Abtasrate erreichen im 
stationären Zustand!? Aber das ist Spekulation von mir.