Определение фазового сдвига между сигналами

Я создал три одинаковые волны с фазовым сдвигом в каждом. Например:

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)

Теперь я хотел бы использовать метод для обнаружения этого сдвига фазы между волнами. Дело в том, что я могу в конечном итоге применить метод к реальным данным и идентифицировать сдвиги фаз между сигналами.

До сих пор я думал о вычислении поперечных спектров между каждой волной и первой волной (т.е. без сдвига фазы):

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

но я не уверен, что это имеет смысл.

Кто-нибудь еще сделал что-то похожее на это? Желаемый результат должен показывать сдвиг фазы на 10 и 15 для волн 2 и 3 соответственно.

Любые советы будут оценены.

Ответ 1

Существует несколько способов измерения фазового сдвига между сигналами. Между вашим ответом, комментариями ниже вашего ответа и другими ответами вы получили большинство вариантов. Конкретный выбор техники обычно основывается на таких проблемах, как:

Шум или чистый: есть ли шум в вашем сигнале?
Многокомпонентный или однокомпонентный: есть ли в вашей записи более одного типа сигнала (несколько тонов на разных частотах, движущихся в разных направлениях)? Или, есть ли только один сигнал, как в примере с синусоидальной волной?
Мгновенный или усредненный: вы ищете среднее фазовое отставание по всей вашей записи, или вы хотите отслеживать, как изменяется фаза во время записи?

В зависимости от вашего ответа на эти вопросы вы можете рассмотреть следующие методы:

Кросс-корреляция: используйте команду a, как [c,lag]=xcorr(y1,y2);, чтобы получить взаимную корреляцию между двумя сигналами. Это работает с исходными сигналами во временной области. Вы ищете индекс, где c является максимальным ([maxC,I]=max(c);), а затем вы получаете свое значение запаздывания в единицах образцов lag = lag(I);. Этот подход дает вам среднее фазовое отставание для всей записи. Это требует, чтобы ваш сигнал, представляющий интерес для записи, был сильнее, чем что-либо еще в вашей записи... другими словами, он чувствителен к шуму и другим помехам.
Частотная область: здесь вы конвертируете свои сигналы в частотную область (используя fft или cpsd или что-то еще). Затем вы найдете бит, который соответствует частоте, о которой вы заботитесь, и получите угол между двумя сигналами. Так, например, если бит # 18 соответствует вашей частоте сигнала, вы получите фазовое отставание в радианах через phase_rad = angle(fft_y1(18)/fft_y2(18));. Если ваши сигналы имеют постоянную частоту, это отличный подход, потому что он, естественно, отвергает все помехи и помехи на других частотах. У вас может быть очень сильное вмешательство на одной частоте, но вы все равно можете получить свой сигнал на другой частоте. Этот метод не самый лучший для сигналов, которые изменяют частоту во время окна анализа fft.
Преобразование Гильберта. Третий метод, который часто упускается из виду, заключается в преобразовании вашего сигнала во временной области в аналитический сигнал через преобразование Гильберта: y1_h = hilbert(y1);. Как только вы это сделаете, ваш сигнал представляет собой вектор комплексных чисел. Вектор, содержащий простую синусоидальную волну во временной области, теперь будет вектором комплексных чисел, величина которых постоянна и фаза меняется синхронно с вашей исходной синусоидальной волной. Этот метод позволяет мгновенно получить фазовый лаг между двумя сигналами... он мощный: phase_rad = angle(y1_h ./ y2_h); или phase_rad = wrap(angle(y1_h) - angle(y2_h));. Основное ограничение этого подхода заключается в том, что ваш сигнал должен быть монокомпонентным, что означает, что ваш сигнал, представляющий интерес, должен доминировать над вашей записью. Поэтому вам может потребоваться отфильтровать любые существенные помехи, которые могут существовать.

Ответ 2

Для двух синусоидальных сигналов фаза сложного коэффициента корреляции дает вам то, что вы хотите. Я могу дать вам пример python (используя scipy), поскольку у меня нет теста Matlab для его тестирования.

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)

В matlab есть функция corrcoeff, которая тоже должна работать (питон отбрасывает мнимую часть). То есть c = corrcoeff (x1h, x2h) должен работать в matlab.

Ответ 3

С правильными сигналами:

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

Следующее должно работать:

>> 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

Является ли это достаточно хорошим для ваших "реальных" сигналов, только вы можете сказать.

Ответ 4

Вот небольшая модификация вашего кода: phi = 10 на самом деле находится в степени, затем в синусоидальной функции информация о фазе в основном выражается в радиан, поэтому вам нужно изменить deg2rad(phi) следующим образом:

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)

затем используя метод частотной области, как упомянуто в фиолете

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)));

теперь это даст вам оценку фазового сдвига с помощью error = +-0.2145

Ответ 5

Код Matlab для поиска относительной фазы с использованием взаимной корреляции:

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