2
$\begingroup$

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)

PLV over time

$\endgroup$
1
  • 1
    $\begingroup$ I'm pretty sure this is a side effect of your hilbert transform. $\endgroup$ – Ben Mar 15 at 16:36
1
$\begingroup$

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

leak

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))
$\endgroup$
12
  • $\begingroup$ I think your code differs from what I did though: As far as I understand, you create a random offset for every timepoint once, and then copy that along the N-times. I created a random offset for every 'trial', but this was then held constant along time. Or did I get something wrong? $\endgroup$ – a_student Mar 16 at 11:45
  • $\begingroup$ @a_student Just to clarify: if your random sequence is [1,5,4], then you have [1, 5, 4; 1, 5, 4; 1, 5, 4] as the array, while I (could) have [1, 5, 4; 2, 0, 3; 2, 9, 8], i.e. only the first row matches, al the others are random. If that's true, then I will update the code, but I still get a single line. $\endgroup$ – a concerned citizen Mar 16 at 13:03
  • $\begingroup$ Sorry if I was not clear before, what I meant is that in my code $x^{k}_{t} = sin(f t + \phi^{k})$, where $k$ is the trial/row, $t$ is the timepoint/column and $\phi^{k}$ is a random phase-offset specific to the trial $k$. Whereas in your (original) code from my understanding $x^{k}_{t} = sin(f t + \phi_{t})$, i.e. a phase-offset specific to the time-point. Anyway, I do not know too much matlab so maybe I am interpreting it wrong. $\endgroup$ – a_student Mar 16 at 14:10
  • $\begingroup$ In any case, it seems strange that there would be such a difference... At least the source code of the python function (github.com/scipy/scipy/blob/v1.6.1/scipy/signal/…) looks exactly like the algorithm described in the Matlab documentation (de.mathworks.com/help/signal/ref/hilbert.html) $\endgroup$ – a_student Mar 16 at 14:11
  • 1
    $\begingroup$ yes, that is what i meant $\endgroup$ – a_student Mar 16 at 17:50

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.