Phasenverschiebung mit Cross Correlation in Matlab

Owen S. OP #3966422
Lesenswert?
• ▲
▼ 
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?

1
hcorr = dsp.Crosscorrelator;

2
 t=[0:0.001:1];

3
 x1=sin(2*pi*2*t);

4
 x2=sin(2*pi*2*t – pi/4);

5
 y=step(hcorr,x2,x1); %computes cross-correlation of x1 and x2

6
 figure,plot(t,x1,'b',t, x2, 'g'); 

7
 legend('Input signal 1',' Input signal 2') 

8
 figure, plot(y); title('Correlated output')


Gruss
Senmeis
Antwort Zitieren Bearbeiten Thread verschieben Thread sperren Anmeldepflicht aktivieren Anpinnen Thread löschen Thread mit anderem zusammenführen Beitrag melden
Frank Gast #3969655
Lesenswert?
• ▲
▼ 
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.
Antwort Zitieren Bearbeiten Löschen Beitrag melden
Owen S. OP #3973232
Lesenswert?
• ▲
▼ 
OK, in diesem Beispiel:
1
t_sample = 0.001;

2
period = 0.5;

3
local_maximum = 157; % Index des ersten Maximums

4
deltaT = t_sample*local_maximum

Die Phase wird wie folgt berechnet:
1
deltaPhi = deltaT*2*pi/period;  % Ergebnis: 1.9729


Das stimmt aber nicht mit pi/4 = 0.7854.

Wo habe ich Fehler gemacht?

Senmeis
Antwort Zitieren Bearbeiten Löschen Beitrag melden
Frank Gast #3973438
Lesenswert?
• ▲
▼ 
Habs mal kurz probiert, sollte einigermaßen passen.

1
clear all;

2
close all;

3


4
%=================

5
%   Einstellungen

6
%=================

7
n=10;               %Anzahl Perioden

8
f=50;               %Signalfrequenz

9
fs=f*1000;         %Abtastfrequenz

10
phiGRAD = 50;       %Gewuenschte Phasenverschiebung in GRAD

11


12
%=================

13
%   Konstanten

14
%=================

15
radToGrad = 180/pi;

16
gradToRad = 1/radToGrad;

17
%=================

18
%   Berechnungen

19
%=================

20
T=1/f;              %Periodendauer

21
ts=1/fs;            %Abtastdauer

22
t=0:ts:(T*n - ts);

23
phiRAD = phiGRAD * gradToRad;

24


25
x1=sin(2*pi*f*t);

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

27


28
corr = conv(x1,x2);             %Faltung

29
tCorr = [fliplr(-t(2:end)) t];  %Zeitvektor fuer Faltung

30


31
h=figure('Name','Kreuzkorreltation zweier Signale');

32
maxCorrAmp = max(abs(corr));                %Absoluter Maximalwert der Kreuzkorrelation

33
[i,j,v] = find(abs(corr) == maxCorrAmp);    %Verschiebung suchen

34
deltaT = tCorr(j+1);                          %Verschiebung berechnen

35


36
phiCalcRad = 2*pi*f*deltaT;                  %aus KK berechnete Phase

37
phiCalcGrad = phiCalcRad * radToGrad;        %in Grad

38
phiErr = abs(phiCalcGrad-phiGRAD);            %Phasen-Fehler

39


40
%=================

41
%   Plots

42
%=================

43
subplot(2,1,1);

44
plot(t,x1,t, x2); 

45
legend('Input signal 1',' Input signal 2');

46
grid on;

47


48
subplot(2,1,2);

49
plot(tCorr,corr);

50
text(tCorr(j), corr(j), strcat('\leftarrow \Deltat = ', num2str(deltaT)), 'FontSize', 10);

51
legend('Crosscorrelation');

52
grid on;

53


54
%=================

55
%   Cmd

56
%=================

57
fprintf('delta_T =  %f s\n',deltaT);

58
fprintf('phiCalcGrad =  %f °\n',phiCalcGrad);

59
fprintf('phiErr =  %f °\n',phiErr);
Antwort Zitieren Bearbeiten Löschen Beitrag melden
Frank Gast #3977031
Lesenswert?
• ▲
▼ 
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...
Antwort Zitieren Bearbeiten Löschen Beitrag melden
Frank Gast #3984224
Lesenswert?
• ▲
▼ 
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.
Antwort Zitieren Bearbeiten Löschen Beitrag melden
Thread beobachten

Antwort schreiben

Bitte melde dich an, um einen Beitrag zu schreiben.

Bestehender Account

oder

Mit Google-Account einloggen

Noch kein Account?

Die Registrierung ist kostenlos und dauert nur eine Minute.

Jetzt registrieren