logo标签列表

识别信号之间的相位差

matlabsignal-processingphase
24

我生成了三个相同的波形，每个波形都有一个相位偏移。例如：

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

我现在想使用一种方法来检测这些波之间的相位差。这样做的目的是为了最终将该方法应用于真实数据，并识别信号之间的相位差。

到目前为止，我一直在考虑计算每个波与第一个波（即没有相位差的波）之间的交叉谱：

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的相移。

任何建议都将不胜感激。

1
请参考http://uk.mathworks.com/matlabcentral/newsreader/view_thread/52460。 - am304
顺便提一下，你代码中的相位移动是10（或15）弧度，而不是10（或15）度。 - am304
也许可以转到信号处理网站？http://dsp.stackexchange.com - Felipe G. Nievinski
7个回答
54

有几种方法可以测量信号之间的相位差。根据您的回答、下面的评论和其他答案，您已经知道了大部分的选择。具体的技术选择通常基于以下问题：

  • 噪声或清晰：信号中是否有干扰噪声？
  • 多组分或单组分：记录中是否存在多种类型的信号 (多个频率上的多个音调在不同方向上移动)？还是仅存在一个单一信号，如正弦波示例中所示？
  • 瞬时或平均：您是要查找整个记录的平均相位滞后，还是要跟踪相位随时间的变化情况？

根据这些问题的答案，您可以考虑以下技术：

  • 交叉相关：使用类似[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));。该方法的主要限制是您的信号必须是单分量的，即您感兴趣的信号必须占据录音中的主导地位。因此，您可能需要滤除任何存在的重要干扰。

    • 非常好的回答。请查看Math.StackExchange上一个问题的我的答案，这是hilbert函数如何用于计算两个信号之间瞬时相对相位的示例。 - horchler
    “lag” 应该是一个函数吗？我找不到它。 - Asad Saeeduddin
    1
    关于“滞后”，它是从“xcorr”输出的。我的原始帖子中有一个打字错误。很抱歉！描述“交叉相关”方法中使用“滞后”的文本现在应该是正确的。 - chipaudette
    @chipaudette FFT分析似乎在相位方面无法正常工作。以下示例使用Python完成：t = np.arange(0,100,step=0.01); f = np.sin(1.5107*np.pi*2*t+np.pi/6); ft = np.fft.fft(f); a=np.where(abs(ft)==np.amax(abs(ft)))[0][1]; A1 = np.rad2deg(np.angle(ft[a])); f = np.sin(1.5254*np.pi*2*t); ft = np.fft.fft(f); a=np.where(abs(ft)==np.amax(abs(ft)))[0][1]; A2 = np.rad2deg(np.angle(ft[a]));预期A2 - A1为30，但实际上是125。只有当两个正弦波具有完全相同的频率时才能正常工作。 - JZYL
    @Jimmy，如果两个信号的频率不完全相同，无论使用哪种分析技术，谈论相位差都有点奇怪。OP假设所有被分析的信号具有相同的频率，因此我也采用了这一点。 - chipaudette
    5

    对于两个正弦信号，复相关系数的相位可以帮助您获得所需结果。由于我没有matlab来测试，我只能提供一个Python示例（使用scipy）。

    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版本会丢弃虚部）。即在matlab中，c = corrcoeff(x1h,x2h)应该可以正常工作。

    2
    我已经搜寻了很久，而这个被埋藏的小答案完美地解决了我的Python脚本问题。不过有几点需要注意！ hilbert 需要 from scipy import signal 才能成为 signal.hilbert；inner、conj 和 angle 需要 import numpy 才能成为 numpy.inner、numpy.conj 和 numpy.angle；sqrt 需要 import math 才能成为 math.sqrt。这对一些人可能很明显，但还是值得澄清的。感谢这个答案！ - erekalper
    如果你稍微扰动正弦频率，它就不再起作用了。例如，将x1的频率更改为0.1010，将x2的频率更改为0.1005。这在实际应用中总是会发生的。我们该如何解决这个问题？ - JZYL
    好的，但如果频率不同，相位就会漂移。那么整个问题原则上就消失了。你唯一能做的就是计算每个窗口的相位偏移量。 - André Bergner
    2

    使用交叉相关寻找相对相位的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
    1

    这是您代码的小修改：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)));

    现在这个代码将会给你一个相位偏移的估计，误差为+-0.2145。
    1
    使用正确的信号：
    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

    无论这是否足够满足你的“真实”信号，只有你自己能够判断。

    0
    0

    如果您使用带延迟的AWGN信号并应用您的方法，它会起作用，但如果您使用单音频率估计将无法帮助您。因为除了该音调外，没有其他频率的能量。您最好在时域中使用互相关来进行处理——对于固定延迟，这样效果更好。如果您有一个宽带信号，可以使用子带域并从中估算相位（由于低交叉频率依赖性，这比FFT更好）。

    网页内容由stack overflow 提供, 点击上面的
    可以查看英文原文，
    原文链接