webentwicklung-frage-antwort-db.com.de

identifizieren der Phasenverschiebung zwischen Signalen

Ich habe drei identische Wellen mit jeweils einer Phasenverschiebung erzeugt. Zum Beispiel:

t = 1:10800; % generate time vector
fs = 1; % sampling frequency (seconds)
A = 2; % amplitude
P = 1000; % period (seconds), the time it takes for the signal to repeat itself
f1 = 1/P; % number of cycles per second (i.e. how often the signal repeats itself every second).
y1 = A*sin(2*pi*f1*t); % signal 1
phi = 10; % phase shift
y2 = A*sin(2*pi*f1*t + phi); % signal 2
phi = 15; % phase shift
y3 = A*sin(2*pi*f1*t + phi); % signal 3

YY = [y1',y2',y3'];

plot(t,YY)

enter image description here

Ich möchte jetzt eine Methode verwenden, um diese Phasenverschiebung zwischen den Wellen zu erkennen. Dies ist so, dass ich die Methode schließlich auf reale Daten anwenden und Phasenverschiebungen zwischen Signalen identifizieren kann. 

Bisher habe ich daran gedacht, die Kreuzspektren zwischen jeder Welle und der ersten Welle zu berechnen (d. H. Ohne Phasenverschiebung):

for i = 1:3;
    [Pxy,Freq] = cpsd(YY(:,1),YY(:,i));
    coP = real(Pxy);
    quadP = imag(Pxy);
    phase(:,i) = atan2(coP,quadP);
end

ich bin mir aber nicht sicher, ob das Sinn macht. 

Hat jemand anderes etwas Ähnliches getan? Das gewünschte Ergebnis sollte eine Phasenverschiebung bei 10 und 15 für die Wellen 2 bzw. 3 zeigen. 

Jeder Rat wäre dankbar. 

17
Emma Tebbs

Es gibt verschiedene Möglichkeiten, die Phasenverschiebung zwischen Signalen zu messen. Zwischen Ihrer Antwort, den Kommentaren unter Ihrer Antwort und den anderen Antworten haben Sie die meisten Optionen erhalten. Die spezifische Wahl der Technik basiert in der Regel auf folgenden Themen: 

  • Laut oder sauber : Gibt es Rauschen in Ihrem Signal? 
  • Multi-Component oder Single-Component : Gibt es mehr als einen Signaltyp in Ihrer Aufnahme (mehrere Töne mit mehreren Frequenzen, die sich in verschiedene Richtungen bewegen)? Oder gibt es nur ein einzelnes Signal, wie in Ihrem Sinuswellenbeispiel? 
  • Momentan oder Durchschnitt : Suchen Sie nach der durchschnittlichen Phasenverzögerung Ihrer gesamten Aufnahme oder möchten Sie verfolgen, wie sich die Phase während der Aufnahme ändert? 

Abhängig von Ihrer Antwort auf diese Fragen können Sie die folgenden Techniken in Betracht ziehen:

  • Kreuzkorrelation : Verwenden Sie den Befehl a wie [c,lag]=xcorr(y1,y2);, um die Kreuzkorrelation zwischen den beiden Signalen zu erhalten. Dies funktioniert mit den ursprünglichen Zeitbereichsignalen. Sie suchen nach dem Index, bei dem c maximal ist ([maxC,I]=max(c);), und dann erhalten Sie Ihren Lag-Wert in Einheiten von Samples lag = lag(I);. Dieser Ansatz gibt Ihnen die durchschnittliche Phasenverzögerung für die gesamte Aufnahme. Dies setzt voraus, dass Ihr interessierendes Signal bei der Aufnahme stärker ist als alles andere bei Ihrer Aufnahme ... mit anderen Worten, es ist empfindlich gegenüber Rauschen und anderen Interferenzen.

  • Frequency Domain : Hier konvertieren Sie Ihre Signale in die Frequenzdomäne (mit fft oder cpsd oder was auch immer). Dann würden Sie den Behälter finden, der der Frequenz entspricht, die Sie interessieren, und den Winkel zwischen den beiden Signalen ermitteln. Wenn zum Beispiel Bin # 18 der Frequenz Ihres Signals entspricht, erhalten Sie die Phasenverzögerung im Radianten über phase_rad = angle(fft_y1(18)/fft_y2(18));. Wenn Ihre Signale eine konstante Frequenz haben, ist dies ein hervorragender Ansatz, da auf natürliche Weise jegliches Rauschen und Interferenzen bei anderen Frequenzen zurückgewiesen wird. Sie können bei einer Frequenz wirklich starke Interferenzen haben, aber Sie können Ihr Signal immer noch sauber auf einer anderen Frequenz erhalten. Diese Technik eignet sich nicht für Signale, deren Frequenz während des FFT-Analysefensters geändert wird.

  • Hilbert-Transformation : Eine dritte oft übersehene Technik besteht darin, Ihr Zeitbereichssignal über die Hilbert-Transformation in ein analytisches Signal umzuwandeln: y1_h = hilbert(y1);. Sobald Sie dies tun, ist Ihr Signal ein Vektor aus komplexen Zahlen. Ein Vektor, der eine einfache Sinuswelle im Zeitbereich enthält, ist jetzt ein Vektor aus komplexen Zahlen, deren Größe konstant ist und deren Phase sich synchron mit Ihrer ursprünglichen Sinuswelle ändert. Diese Technik ermöglicht es Ihnen, die augenblickliche Phasenverzögerung zwischen zwei Signalen zu erhalten ... es ist mächtig: phase_rad = angle(y1_h ./ y2_h); oder phase_rad = wrap(angle(y1_h) - angle(y2_h));. Die größte Einschränkung bei diesem Ansatz besteht darin, dass Ihr Signal einkomponentig sein muss, was bedeutet, dass Ihr interessierendes Signal Ihre Aufnahme dominieren muss. Daher müssen Sie möglicherweise alle möglicherweise vorhandenen Interferenzen herausfiltern.

44
chipaudette

Bei zwei sinusförmigen Signalen gibt die Phase des komplexen Korrelationskoeffizienten an, was Sie wollen. Ich kann Ihnen nur ein Python-Beispiel geben (mit scipy), da ich kein Matlab habe, um es zu testen.

x1 = sin( 0.1*arange(1024) )
x2 = sin( 0.1*arange(1024) + 0.456)
x1h = hilbert(x1)
x2h = hilbert(x2)
c = inner( x1h, conj(x2h) ) / sqrt( inner(x1h,conj(x1h)) * inner(x2h,conj(x2h)) )
phase_diff = angle(c)

Es gibt eine Funktion corrcoeff in matlab, die auch funktionieren sollte (der Python man den imaginären Teil verwerfen). D.h. c = corrcoeff (x1h, x2h) sollte in Matlab funktionieren.

3
André Bergner

Der Matlab-Code zum Finden der relativen Phase unter Verwendung der Kreuzkorrelation:

fr = 20; % input signal freq
timeStep = 1e-4;
t = 0:timeStep:50; % time vector
y1 = sin(2*pi*t); % reference signal
ph = 0.5; % phase difference to be detected in radians
y2 = 0.9 * sin(2*pi*t + ph); % signal, the phase of which, is to be measured relative to the reference signal

[c,lag]=xcorr(y1,y2); % calc. cross-corel-n
[maxC,I]=max(c); % find max
PH = (lag(I) * timeStep) * 2 * pi; % calculated phase in radians

>> PH

PH =

    0.4995
1
user3804598

Mit den richtigen Signalen:

t = 1:10800; % generate time vector
fs = 1; % sampling frequency (seconds)
A = 2; % amplitude
P = 1000; % period (seconds), the time it takes for the signal to repeat itself
f1 = 1/P; % number of cycles per second (i.e. how often the signal repeats itself every second).
y1 = A*sin(2*pi*f1*t); % signal 1
phi = 10*pi/180; % phase shift in radians
y2 = A*sin(2*pi*f1*t + phi); % signal 2
phi = 15*pi/180; % phase shift in radians
y3 = A*sin(2*pi*f1*t + phi); % signal 3

Folgendes sollte funktionieren:

>> acos(dot(y1,y2)/(norm(y1)*norm(y2)))
>> ans*180/pi
ans =  9.9332
>> acos(dot(y1,y3)/(norm(y1)*norm(y3)))
ans =  0.25980
>> ans*180/pi
ans =  14.885

Ob das für Ihre "echten" Signale gut genug ist oder nicht, können nur Sie sagen.

1
am304

Hier ist die kleine Modifikation Ihres Codes: phi = 10 ist tatsächlich in Grad, dann werden die Phaseninformationen in der Sinus-Funktion meistens im Bogenmaß ausgedrückt. Sie müssen deg2rad(phi) wie folgt ändern:

t = 1:10800; % generate time vector
fs = 1; % sampling frequency (seconds)
A = 2; % amplitude
P = 1000; % period (seconds), the time it takes for the signal to repeat itself
f1 = 1/P; % number of cycles per second (i.e. how often the signal repeats itself every second).
y1 = A*sin(2*pi*f1*t); % signal 1
phi = deg2rad(10); % phase shift 
y2 = A*sin(2*pi*f1*t + phi); % signal 2
phi = deg2rad(15); % phase shift
y3 = A*sin(2*pi*f1*t + phi); % signal 3

YY = [y1',y2',y3'];

plot(t,YY)

dann unter Verwendung der Frequenzbereichsmethode, wie Chipaudette erwähnt

fft_y1 = fft(y1);
fft_y2 = fft(y2);
phase_rad = angle(fft_y1(1:end/2)/fft_y2(1:end/2));
phase_deg = rad2deg(angle(fft_y1(1:end/2)/fft_y2(1:end/2)));

jetzt erhalten Sie eine Phasenverschiebung mit error = +-0.2145.

1
pcu