Выявление фазового сдвига между сигналами

Question

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

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 соответственно.

Любой совет будет принят во внимание.

Answer 1

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

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

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

• Кросс-корреляция: используйте команду как [c,lag]=xcorr(y1,y2); чтобы получить взаимную корреляцию между двумя сигналами. Это работает с исходными сигналами во временной области. Вы ищите индекс, где c максимум ([maxC,I]=max(c);), а затем вы получите значение лага в единицах выборок lag = lag(I);, Этот подход дает вам среднее запаздывание фазы для всей записи. Это требует, чтобы ваш интересующий вас сигнал в записи был сильнее, чем что-либо еще в вашей записи... другими словами, он чувствителен к шуму и другим помехам.

• Частотная область: здесь вы преобразуете свои сигналы в частотную область (используя fft или же cpsd или что угодно). Затем вы найдете корзину, соответствующую частоте, которая вас интересует, и получите угол между двумя сигналами. Так, например, если bin #18 соответствует частоте вашего сигнала, вы получите фазовый лаг в радианах через phase_rad = angle(fft_y1(18)/fft_y2(18));, Если ваши сигналы имеют постоянную частоту, это отличный подход, поскольку он естественным образом подавляет все помехи и помехи на других частотах. Вы можете иметь очень сильные помехи на одной частоте, но вы все равно можете получить чистый сигнал на другой частоте. Этот метод не является лучшим для сигналов, которые изменяют частоту во время окна анализа БПФ.

• Преобразование Гильберта. Третий метод, который часто упускают из виду, заключается в преобразовании вашего сигнала во временной области в аналитический сигнал с помощью преобразования Гильберта: y1_h = hilbert(y1);, Как только вы это сделаете, ваш сигнал будет вектором комплексных чисел. Вектор, содержащий простую синусоидальную волну во временной области, теперь будет вектором комплексных чисел, величина которого постоянна и чья фаза изменяется синхронно с вашей исходной синусоидой. Эта техника позволяет получить мгновенный фазовый интервал между двумя сигналами... он мощный: phase_rad = angle(y1_h ./ y2_h); или же phase_rad = wrap(angle(y1_h) - angle(y2_h));, Основным ограничением этого подхода является то, что ваш сигнал должен быть однокомпонентным, а это означает, что ваш интересующий сигнал должен доминировать в вашей записи. Поэтому вам, возможно, придется отфильтровать любые существенные помехи, которые могут существовать.

Answer 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, которая тоже должна работать (Python отбрасывает мнимую часть). Т.е. c = corrcoeff (x1h, x2h) должен работать в Matlab.

Answer 3

Код 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

Answer 4

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

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

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

Answer 5

Вот небольшая модификация вашего кода: 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

Answer 6

Если вы используете сигнал AWGN с задержкой и применяете свой метод, он работает, но если вы используете однотональную оценку частоты, вам это не поможет. потому что на любой другой частоте, кроме тона, нет энергии. Для этого лучше использовать взаимную корреляцию во временной области - она ​​лучше работает при фиксированной задержке. Если у вас есть широкополосный сигнал, вы можете использовать область поддиапазонов и оценивать фазу по нему (это лучше, чем БПФ из-за низких межчастотных зависимостей).

Answer 7

Если вы знаете частоту и просто хотите найти фазу, а не использовать полное БПФ, вы можете рассмотреть алгоритм Герцеля, который является более эффективным способом вычисления ДПФ для одной частоты (БПФ вычислит его для все частоты).

Для хорошей реализации см: https://www.mathworks.com/matlabcentral/fileexchange/35103-generalized-goertzel-algorithm и https://asp-eurasipjournals.springeropen.com/track/pdf/10.1186/1687-6180-2012-56.pdf