为什么锁相值会在信号的开头和结尾显示杂散效应？

信息处理 信号分析 阶段 时域
2022-02-09 12:38:19

我试图通过一些简单的例子来理解锁相值的行为，例如定义在这里。

基本上，我正在创建两个具有相同频率和幅度的信号，其中两者之间的相位偏移是随机的。现在，正如预期的那样，在大多数时间点，PLV 给了我一个低值，表明没有一致的相位开始模式。然而，在信号的开头和结尾，有一些奇怪的效果。我想这些可能是由希尔伯特变换以某种方式引入的，还是我的实现有问题？

这是我定义函数的方式：

import numpy as np
from scipy import signal
import matplotlib.pyplot as plt


def plv(x,y,trial_axis=0):

    hilbert_x = signal.hilbert(x)
    hilbert_y = signal.hilbert(y)
    
    normed_x = hilbert_x/np.absolute(hilbert_x)
    normed_y = hilbert_y/np.absolute(hilbert_y)
    
    
    p = np.absolute(np.mean(normed_x*normed_y.conj(),axis=trial_axis))
    
    return p

在这里，我们有一个行为图：

t = np.linspace(0,10,1000)
N = 50
s_1 = np.array([np.sin(10*2*np.pi*t) for i in range(N)])
s_2 = np.array([np.sin(10*(2*np.pi*t)+np.random.uniform(low=0,high=2*np.pi) ) for i in range(N)])
p = plv(s_1,s_2)
plt.plot(p)

PLV 随时间变化

1个回答

以前作为答案发布的内容应该是扩展评论——对此感到抱歉。根本没有足够的空间来写我想写的关于 Octave 代码的行为方式的内容。既然代码已被清除，那么首先应该是答案（鉴于相当明显的迹象），现在答案是：频谱泄漏。

希尔伯特变换是通过将信号的一半 FFT 归零然后通过 IFFT 重新计算分析信号来计算的。由于数据长度有限的性质，蛮力归零总是会导致一些泄漏。这可以很容易地用您的代码进行测试。s1是一个简单的sin()，所以它的希尔伯特变换将是一个-cos()。只需绘制 之间的差异-cos()-hx.imag，您就会看到泄漏。第二个信号也是如此。

这是 Python 中更改后的代码：

import numpy as np
from scipy import signal
import matplotlib.pyplot as plt
f = 10
t = np.linspace(0, 100/f, 1000)
N = 50f = 10
nt = 1000
t = np.linspace(0, 100/f, nt)
N = 50
hor = np.ones(nt)
ver = np.ones(N)
r = np.array( [ np.random.uniform(low=0, high=2*np.pi) for i in range(N) ] )
s1 = np.array( [ np.sin( 2*np.pi*f*t ) for i in range(N) ] )
s2 = np.array( np.sin( 2*np.pi*( f*ver[:,None]*t + hor*r[:,None] ) ) )
c1 = -np.array( [ np.cos(2*np.pi*f*t) for i in range(N) ] )
c2 = -np.array( np.cos( 2*np.pi*( f*ver[:,None]*t + hor*r[:,None] ) ) )
hx = signal.hilbert(s1)
hy = signal.hilbert(s2)
nx = hx/np.absolute(hx)
ny = hy/np.absolute(hy)
p = np.absolute( np.mean( nx*ny.conj(), axis=0 ) )
plt.plot(p)
plt.show()

这是测试（c2由于随机性将显示非重叠图）：

plt.plot(c1.transpose() - hx.imag.transpose())
plt.show()

泄漏

Octave 代码，用于检查：

f = 10
t = linspace(0, 100/f, 1000);
N = 50;
c = ones(N, 1);
ft = f*t.*c;
s1 = sin(2*pi*ft);
c1 = -cos(2*pi*ft);
r = rand(N,1).*c;
s2 = sin(2*pi*(ft + r));
c2 = -cos(2*pi*(ft + r));
hx = hilbert(s1');
hy = hilbert(s2');
nx = hx./abs(hx);
ny = hy./abs(hy);
p = abs(mean((nx.*conj(ny))'));
grid on
subplot(2, 1, 1)
plot(p)
subplot(2, 1, 2)
plot(c1' - imag(hx))
其它你可能感兴趣的问题
上一篇如何量化两个信号的频谱信息差异？ 下一篇卡尔曼滤波器 - 导出状态转换函数