识别信号之间的相移答案

【问题标题】：identifying phase shift between signals识别信号之间的相移
【发布时间】：2015-02-17 03:46:20
【问题描述】：

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

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

但我不确定这是否有意义。

有没有其他人做过类似的事情？预期的结果应该分别显示第 2 波和第 3 波在 10 和 15 处的相移。

任何建议将不胜感激。

【问题讨论】：

参见例如uk.mathworks.com/matlabcentral/newsreader/view_thread/52460
顺便说一句，您的代码中的相移是 10 (resp. 15) rad，而不是 10 (resp. 15) deg
可能会进入信号处理站点？ dsp.stackexchange.com

标签： matlab signal-processing phase

【解决方案1】：

您可以通过多种方式测量信号之间的相移。在您的回复、您的回复下方的 cmets 以及其他答案之间，您已经获得了大部分选项。技术的具体选择通常基于以下问题：

嘈杂或干净：您的信号中有噪音吗？
多分量或单分量：您的录音中是否存在不止一种类型的信号（多个频率的多个音调向不同方向移动）？或者，是否只有一个信号，例如您的正弦波示例？
瞬时或平均：您是在寻找整个录制过程中的平均相位滞后，还是想要跟踪整个录制过程中的相位变化？

根据您对这些问题的回答，您可以考虑以下技巧：

Cross-Correlation：使用[c,lag]=xcorr(y1,y2); 之类的命令来获取两个信号之间的互相关。这适用于原始时域信号。您查找c 最大的索引（[maxC,I]=max(c);），然后以样本lag = lag(I); 为单位获得滞后值。这种方法为您提供了整个记录的平均相位滞后。它要求您在录音中感兴趣的信号比录音中的任何其他信号都强……换句话说，它对噪音和其他干扰很敏感。
频域：在这里您将信号转换为频域（使用fft 或cpsd 或其他）。然后，您会找到与您关心的频率相对应的 bin，并获得两个信号之间的角度。因此，例如，如果 bin #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));。这种方法的主要限制是您的信号必须是单分量的，这意味着您感兴趣的信号必须在您的录音中占主导地位。因此，您可能必须过滤掉可能存在的任何实质性干扰。

【讨论】：

很好的答案。有关如何使用 hilbert 函数计算两个信号之间的瞬时相对相位的示例，请参阅 Math.StackExchange 上的问题 my answer。
lag 应该是一个函数吗？我在任何地方都找不到它。
关于lag，是xcorr的输出。我原来的帖子有错别字。对不起！描述在“互相关”方法中使用lag 的文本现在应该是正确的。
@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。只有当两个正弦曲线具有精确相同的频率时才有效。
@Jimmy，如果两个信号的频率不完全相同，那么无论采用何种分析技术，谈论相位差都会有点奇怪。 OP 假设所有正在分析的信号的频率相同，所以我同意了。

【解决方案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 丢弃了虚部）。 IE。 c = corrcoeff(x1h,x2h) 应该可以在 matlab 中使用。

【讨论】：

我一直在寻找很长时间，这个隐藏的小答案非常适合我的 Python 脚本。不过要注意几点！ hilbert 需要 from scipy import signal 才能实现 signal.hilbert; inner、conj 和 angle 需要 import numpy 才能使它们成为 numpy.inner、numpy.conj 和 numpy.angle； sqrt 需要 import math 才能使其成为 math.sqrt。这些对某些人来说可能是显而易见的，但需要澄清。感谢您的回答！
如果您稍微干扰正弦频率，这将不再有效。例如，将 x1 的频率更改为 0.1010，将 x2 的频率更改为 0.1005。对于实际应用程序，情况总是如此。我们将如何解决这个问题？
好的，但是如果频率不同，那么相位会漂移。然后整个问题原则上消失了。你唯一能做的就是计算每个窗口的相移。

【解决方案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

【讨论】：

【解决方案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

这对于您的“真实”信号是否足够好，只有您自己知道。

【讨论】：

【解决方案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)

然后使用提到的频域方法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 的相移估计

【讨论】：

【解决方案6】：

如果您知道频率并且只想找到相位，而不是使用完整的 FFT，您可能需要考虑 Goertzel 算法，这是一种计算单个频率的 DFT 的更有效方法（FFT 将计算所有频率）。

如需良好的实现，请参阅：https://www.mathworks.com/matlabcentral/fileexchange/35103-generalized-goertzel-algorithm 和 https://asp-eurasipjournals.springeropen.com/track/pdf/10.1186/1687-6180-2012-56.pdf

【讨论】：

【解决方案7】：

如果您使用带延迟的 AWGN 信号并应用您的方法，它会起作用，但如果您使用单音频率估计将无济于事。因为除了音调之外，任何其他频率都没有能量。为此，您最好在时域中使用互相关 - 它对于固定延迟会更好。如果您有宽带信号，则可以使用子带域并从中估计相位（由于交叉频率依赖性低，它比 FFT 更好）。

【讨论】：