非线性系统 - 计算两个正弦波之间的相位差 - 傅里叶和希尔伯特

Question

编辑：从建议使用xcorr函数得出的结果为零来看，这些数据更可能代表非线性系统，而不是我最初想的那样。在这种情况下，a和b的表达式被简化了，不一定完全与真实数据相比。此外，相移d的值不能再被假定为常数。

我有两个原始数据数组，对应于同一（已知）频率上的两个正弦信号随时间变化，其中一个信号相位被未知量delta，d所改变。这些信号可以描述为：

a = A*sin(wt)
b = B*sin(wt+d)

其中w = 2*pi*f，A, B分别为a和b的振幅。

这里的目标是评估给定信号随时间的相位移动。这可以作为瞬时值或在许多周期内平均值来完成（相对于总周期数很小，即在f = 150Hz下，10秒测试相当于1500个周期）。

我知道有许多方法可以用来解决这个问题，许多其他用户已经就此提出/回答了问题。下面是链接到原始帖子和我实现的代码的方法。我目前遇到的问题是：

1. 如何解释评估后的相移（对于我的情况，相移通常被引用为一个单一的、正数的数字，在许多周期内保持相对稳定）
2. 不同方法之间的值差异
3. 是否有任何我错过的额外方法？

对于我的数据，a和b都是[1 x 399]的大小，涵盖了20个周期。这只是预期数据的一小部分，我正在使用它来测试不同的方法。可以合理地期望至少有10度的相位移，尽管确切数值未知：

傅里叶变换

来源: 如何计算两个正弦信号之间的相位差？

fft_a = fft(a);
fft_a = fft_a(2:end);       %Drop the first point
angle_a = angle(fft_a);     %Angle the result

fft_b = fft(b);
fft_b = fft_b(2:end);       %Drop the first point
angle_b = angle(fft_b);     %Angle the result

ps1 = rad2deg(abs(angle_a - angle_b));     %Phase shift calculation

使用这种方法，我该如何消除开头/结尾的峰值？我需要设置变换在一定数量的点上执行吗？

然后，我该如何绘制相移与时间的图像，因为这个结果现在肯定是在频域中？

希尔伯特变换

来源：识别信号之间的相移

ha = hilbert(a);                 %Hilbert transform
hb = hilbert(b);

ps2 = rad2deg(angle(hb./ha));    %Phase shift calculation

这里得到的相位移随时间变化既有正值又有负值，并呈波形振荡。它的形状相当一致，但我不知道如何解释结果。在对相位移进行 abs 值处理后，结果现在在 0 到 30 度之间变化。

图表： 对于这两种方法，我该如何使用结果来引用在一组周期内的单个相移值？取 mean 看起来不是一个精确的方法。

快速计算相位移我使用了：

ps3 = rad2deg(acos(dot(a,b)/(norm(a)*norm(b))))    %Norm product method

ps3 =

    11.8289

我对这种分析方法还比较新，并且对Matlab不太熟悉，所以非常感谢您提供的任何更正/指导。

编辑：鉴于这是一个非线性系统，有哪些方法可以最好地评估相位移？

- Natalie

---

你看过fminsearch吗？只需使用可参数化的相位偏移减去一个信号，优化最小平方差之和。 - Rody Oldenhuis

...但是对于噪声比较敏感。此外，如果信号振幅不同，需要先进行归一化处理。 - Rody Oldenhuis

这是使用xcorr函数的替代方法吗？此外，似乎我可以对每对点评估相位偏移，而不是整个数组的相位偏移......感谢@rodyoldenhuis的建议，我一定会尝试它。 - Natalie

Answer 1

交叉相关如何处理？

使用交叉相关，您应该能够找到一个绝对的最大偏移值，因为当您将其中一个信号正确地移动时，两个信号将变得相同。

a = A*sin(wt)
b = B*sin(wt+d)
[c, time_c] = xcorr(a,b)
[m, i_m] = max(c)
d_ext = time_c(i_m) * T * w % T is the sampling time

在常量 d 的情况下应该可以工作，但对于变量情况我不确定。

- bracco23

---

我刚刚尝试使用xcorr函数，它给了我一个零的答案......这是否意味着相位移存在，但通过使用此函数被平均为零？我对这个函数或它的操作不太熟悉，无法确定为什么会得到零的答案。 - Natalie

该函数仅实现了交叉相关，您可以在wikipedia页面上找到更多信息。什么是零？您尝试过上面的代码吗？当我尝试使用一个恒定的偏移量时，它能够正常工作。 - bracco23

我认为这更多与我使用的数据类型有关，而不是哪种方法“更好”...是的，你说得对，如果d是一个常量值，它应该可以正常工作。 - Natalie

在您所提供的维基百科页面上，有一个关于非线性系统和零值结果的有趣段落：“……这个问题的出现是因为某些二次矩可以等于零，这可能会错误地表明两个信号之间几乎没有“相关性”（在统计依赖性的意义上），而实际上这两个信号通过非线性动力学密切相关。” 因此，可以合理地得出结论，这里存在非线性相移。对此有何想法？ - Natalie

需要注意的是，xcorr 是 信号处理工具箱 的一部分。 - Rody Oldenhuis

Answer 2

我的意思是：

clc

% Noiseless example data
x = 0:0.1:40*pi;

p1 = deg2rad(+5.74);
p2 = deg2rad(-4.28);

y = sin(x + p1);
z = sin(x + p2);

input = rad2deg(p1 - p2) % display injected phase difference


% Retrieve values using no assumptions 

options = optimset('TolX'   , 0,...
                   'TolFun' , 0, ...
                   'display', 'off');

p1 = fminsearch(@(q) sum( (y - sin(x + q)).^2 ), 0, options);
p2 = fminsearch(@(q) sum( (z - sin(x + q)).^2 ), 0, options);

computed = rad2deg(p1 - p2) % display computed phase difference

在此示例中，输出结果为：

input =
    1.002000000000000e+001
computed =
    1.002000000000000e+001

但是，正如所指出的那样，这种方法相当简单，当噪声出现时，精度下降得比其他方法快得多。

- Rody Oldenhuis