matlab - 识别信号之间的相移

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

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

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

任何意见，将不胜感激。

Answer 1

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

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

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

• 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));. 这种方法的主要限制是您的信号必须是单分量的，这意味着您感兴趣的信号必须在您的录音中占主导地位。因此，您可能必须过滤掉可能存在的任何实质性干扰。

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)

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

Answer 6

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

Answer 7

如果您知道频率并且只想找到相位，而不是使用完整的 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