идентификация сдвига фаз между сигналами

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

Любой совет был бы признателен.

878 5

signal-processing matlab phase

5 ответов:

Comments

Ничего не найдено.

chipaudette · Accepted Answer · 2015-04-06 20:53:40

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

шумно или Чисто : есть ли шум в вашем сигнале?

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

мгновенный или усредненный: вы ищете среднюю фазовую задержку по всей записи, или вы хотите отслеживать, как изменяется фаза на протяжении всей записи?

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

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

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

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

André Bergner · Accepted Answer · 2014-12-19 11:22:57

Для двух синусоидальных сигналов фаза комплексного коэффициента корреляции дает вам то, что вы хотите. Я могу только привести пример 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 One отбрасывает мнимую часть). То есть c = corrcoeff(x1h,x2h) должен работать в matlab.

am304 · Accepted Answer · 2014-12-18 13:42:06

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

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

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

pcu · Accepted Answer · 2016-01-14 04:56:30

Вот небольшая модификация вашего кода: 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)
Затем, используя метод частотной области, как упоминалось chipaudette
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

user3804598 · Accepted Answer · 2017-12-14 14:59:11

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