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

Kreuzkorrelation zweier Signale

 

Thomas K
Forum-Anfänger

Forum-Anfänger


Beiträge: 47
Anmeldedatum: 15.04.12
Wohnort: ---
Version: ---
     Beitrag Verfasst am: 24.07.2013, 15:15     Titel: Kreuzkorrelation zweier Signale
  Antworten mit Zitat      
Hallo,

ich muss die Verschiebung zweier Signale herausfinden, und da habe ich zuerst einmal an eine Kreuzkorrelation der beiden Signale gedacht. Da das für mich Neuland ist, habe ich versucht mich an einem einfachen Testbeispiel einzuarbeiten, nämlich die Phasenverschiebung zwischen zwei Sinus-Signalen zu finden:
Code:

phi = linspace(0,2*pi,360);

x = sin(phi)+1;
y = sin((phi-pi/3));
figure;plot(phi,x,'r');hold all;plot(phi,y,'b');

[val lag] = xcorr(x,y,'coeff');

figure;plot(lag,abs(val));

[m i] = max(val);

theta = phi(abs(lag(i)));

degree = theta*180/pi;
 


Trotz der sehr einfachen Signale komme ich leider damit nur in die Nähe der tatsächlichen, von mir vorgegebenen Phasenverschiebung. Die Abweichung beträgt zwischen 5 und 15 Grad, dass ist einfach nicht brauchbar (1 Grad wäre mir recht).

Da ich nicht weiß, wie Matlab diese Kreuzkorrelation durchführt, habe ich zum Vergleich ein eigenes kleines Skript dazu verfasst (auf einfachste Art und Weise):

Code:

close all;
clear var;
phi = linspace(0,2*pi,360);

x = sin(phi);
y = sin((phi-pi/6));
figure;plot(phi,x,'r');hold all;plot(phi,y,'b');

corr = zeros(size(x,1),1);
sXX = sum(x.^2);

for k=1:size(x,1)
    sXY = 0;
    sYY = 0;
    for i=1:size(x,1)-k
        sXY = sXY+x(i)*y(i+k);
        sYY = sYY+y(i+k)^2;
    end

    corr(k) = sXY/sqrt(sXX*sYY);
end

[maxi ind] = max(corr);
theta = phi(ind);
 


Trotz der paar primitiven Zeilen, sind die Ergebnisse der selbsterstellten Version besser als jene mit Matlab-Funktion xcorr() (ist aber wesentlich langsamer als xcorr()). Hat jemand Vorschläge, wo man da ansetzen könnte?

Schöne Grüße, Thomas.
Private Nachricht senden Benutzer-Profile anzeigen


Winkow
Moderator

Moderator



Beiträge: 3.842
Anmeldedatum: 04.11.11
Wohnort: Dresden
Version: R2014a 2015a
     Beitrag Verfasst am: 24.07.2013, 15:54     Titel:
  Antworten mit Zitat      
hallo thomas.
ich muss leider gestehen das ich sowas noch nicht gemacht hab aber ich hab 2 links gefunden die auf den ersten blick so aussehen als könnten sie dir helfen
http://homepages.hs-bremen.de/~krausd/iwss/SA_Feng_Helber.pdf
und hier im forum http://www.gomatlab.de/frequenz-pha.....-zwei-signalen-t5388.html der beitrag von Maddy.
ich hoffe jemand anders kann dir noch direkt weiterhelfen
grüße winkow
_________________

richtig Fragen
Private Nachricht senden Benutzer-Profile anzeigen
 
Thomas K
Themenstarter

Forum-Anfänger

Forum-Anfänger


Beiträge: 47
Anmeldedatum: 15.04.12
Wohnort: ---
Version: ---
     Beitrag Verfasst am: 24.07.2013, 15:57     Titel:
  Antworten mit Zitat      
Danke für die schnelle Antwort. Werde mir diese Links gleich mal anschauen.
Private Nachricht senden Benutzer-Profile anzeigen
 
bastian79
Forum-Anfänger

Forum-Anfänger


Beiträge: 13
Anmeldedatum: 17.07.13
Wohnort: ---
Version: ---
     Beitrag Verfasst am: 25.07.2013, 09:05     Titel:
  Antworten mit Zitat      
Wenn man wissen will, was innerhalb einer Funktion passiert .... dann mach sie selber!

Code:

phi = linspace(0,2*pi,360);

x = sin(phi)+1;
y = sin((phi-pi/3));

X = fft(x);
Y = fft(y);
C  = conj(X) .* Y;

PHASE = angle(C);
 


gruß
b
Private Nachricht senden Benutzer-Profile anzeigen
 
Thomas K
Themenstarter

Forum-Anfänger

Forum-Anfänger


Beiträge: 47
Anmeldedatum: 15.04.12
Wohnort: ---
Version: ---
     Beitrag Verfasst am: 06.08.2013, 09:24     Titel:
  Antworten mit Zitat      
Hallo noch einmal,

habe mich leider längere Zeit nicht damit befassen können. Ich habe probiert die Zeitverschiebung zweier einfacher Test- Signale (basierend auf dem Beispiel der Matlab-Hilfe zu fft() ) zu detektieren. Dazu nutzte ich einen Algorithmus, den ich in einem von Winkows Pdfs gefunden habe, hier das Pdf:

http://homepages.hs-bremen.de/~krausd/iwss/SA_Feng_Helber.pdf

Ich habe versucht, den Algorithmus umzusetzen, welcher die Fourier-Transformation nutzt, um die Verschiebung im Zeitbereich zu berechnen. Allerdings bekomme ich nicht die richtigen Werte, nicht einmal die Richtung der Verschiebung stimmt. Ich poste hier mal den Code:

Code:

function sol = timeDelay()
close all;
%
% Generate two signals y1 and y2:
%
Fs = 1024;                    %Sampling frequency
T = 1/Fs;                     %Sample time
L = 1024;                     %Length of signal
t1 = (0:L-1)*T;               %Time vector
%
% Signal 1:
%
y1 = 0.7*sin(2*pi*50*t1) + sin(2*pi*120*t1);    %Sum of a 50 Hz sinusoid and a 120 Hz sinusoid
% y1 = x1 + 2*randn(size(t1));     %Sinusoids plus noise
figure;
plot(Fs*t1(1:50),y1(1:50))
title('Signal Corrupted with Zero-Mean Random Noise')
xlabel('time (milliseconds)')
%
delay = 1/1024*11    %Time delay between signal 1 and 2
%
% Signal 2:
%
y2 = 0.7*sin(2*pi*50*(t1+delay)) + sin(2*pi*120*(t1+delay));
% y2 = x2 + 2*randn(size(t1));     %Sinusoids plus noise
figure;
plot(Fs*t1(1:50),y2(1:50))
title('Signal Corrupted with Zero-Mean Random Noise')
xlabel('time (milliseconds)')
%
% Calculate time delay:
%
tP = 1/1024;                %Sampling period
%
m = 2^3;                    %number of data blocks
n = length(y1)/m;           %number of values in each data block
%
Y1 = fft(reshape(y1,n,m));  %FFT of y1
Y2 = fft(reshape(y2,n,m));  %FFT of y2
%
Y1k = conj(Y1);             %Conjugated complex FFT
Y2k = conj(Y2);             %Conjugated complex FFT
%
C11 = mean(Y1.*Y1k,2);      %Auto correlation Y1k
C22 = mean(Y2.*Y2k,2);      %Auto correlation Y2k
C12 = mean(Y1.*Y2k,2);      %Cross correlation Y1k and Y2k
%
Q12 = abs(C12./sqrt(C11.*C22)).^2;  %Square of coherence function
%
W12 = Q12./(1-Q12);         %Weight function
%
% calculate start value for fminsearch:
%
u = unwrap(angle(C12(1:n/4)));
v = (0:n/4-1)'*2*pi/(n*tP);
tStart = -(n/4*sum(v.*u)-sum(v)*sum(u))/(n/4*sum(v.*v)-sum(v)^2);
%
ts = fminsearch(@(ts) opti(ts,W12(1:n/4),unwrap(angle(C12(1:n/4)))),-2*pi*tStart/(n*tP));
%
% Estimation of time delay:
%
sol = -ts*n*tP/(2*pi);           %time delay
%
dT = subTimeDelay(y1,y2,T)        %time delay by subfunction subTimeDelay
%
end
%
% Subfunction opti:
%
function y = opti(ts,W,angleC12)
y = -sum(W.*cos(angleC12-(ts*(0:length(W)-1)')));
end
%
% Subfunction:
function dT = subTimeDelay(y1,y2,T)
% [c l] = xcorr(y1,y2,'coeff');   %Cross correlation (normalized)
% figure;plot(l,c);               %Plot coefficients and lags
l = length(y1);
nFFT = 2^nextpow2(l);       %Next power of 2 from length of y
Y1 = fft(y1,nFFT);          %Fast fourier transform of signal 1
Y2 = fft(y2,nFFT);          %Fast fourier transform of delayed signal 2
dT = zeros(length(Y1),1);

k = 0:1:l/2-1;
for i=1:length(k)
    dT(i) = 1024*T/(2*pi*k(i))*angle(Y1(i)*conj(Y2(i)));    %Time delay with respect to equation (3.3.29) from Pdf
end

ind = find(~isnan(dT));
dT = dT(ind);

end
 


Die letzte Subfunction ist eine sehr reduzierte Version. Da habe ich versucht über den Verschiebungssatz möglichst einfach die Zeitverschiebung zu berechnen. Hat aber auch nicht geklappt.

Vielleicht hat jemand Erfahrung damit, und weiß was anders gemacht werden muss. Ich wird mich damit auch noch etwas beschäftigen.

Schöne Grüße, Thomas.
Private Nachricht senden Benutzer-Profile anzeigen
 
Neues Thema eröffnen Neue Antwort erstellen



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.