Why does Phase-locking value show spurious effects at the beginning & end of a signal?

Question

I am trying to understand the behaviour of the phase-locking-value, as e.g. defined here, by some simple examples.

Basically, I am creating two signals with the same frequency & amplitude, where the phase-offset between the two is random. Now, as expected, at most time points the PLV gives me a low value indicating that there is not a consistent pattern of phase onset. However, at the beginning & end of the signal, there are some strange effects. I suppose that these might be somehow introduced by the Hilbert transform, or is this a problem with my implementation?

Here is how I defined the function:

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

Here we have a plot of the behaviour:

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)

Answer 1

What was previously posted as an answer should have been an extended comment -- sorry for that. There simply was not enough space for what I thought I'd write about how the Octave code seemed to behave. Now that the code is cleared, what should have been the answer in the first place (given the rather obvious signs), is now the answer: spectral leakage.

The Hilbert transform is calculated by zeroing half the FFT of the signal and then recomputing the analytical signal through the IFFT. Brute-force zeroing will always lead to some leakage due to the nature of the limited length of the data. This can be easily tested with your code. s1 is a simple sin(), so its Hilbert transform will be a -cos(). Simply plot the difference between -cos()-hx.imag, and you'll see the leakage. The same goes for the second signal.

This is the changed code in 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()

and this is the test (c2 will show non-overlapping plots due to the randomness):

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

The Octave code, for checking:

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))