LFP-LFP pairwise phase consistency (PPC) is an extension of idea of previously defined spike-field PPC analysis. Here, each of the 20 theta phases in each individual theta cycle was treated as an 'spike event' similar as in spike-field PPC. The following Step1-3 express main idea and details for implementation. One can follow the main idea and implement yourself. In this paper, I used wavelet method and defualt functions from Mathwork, which depends on a lot on matlab versions and differs with wavelet type. One can also use hilber transform of filtered LFP to get phase for a specified frequency band.

Step1: calculate wavelet cross-spectrum (WCS) between two LFPs.

Noted that the cross-spectrum is complex number, whose angle is the phase. The result would be a matrix in frequency * time domain, in the paper there is 81 frequeicies points in frequency domain (20:2:180 Hz).

Step2: Sepearte cross-spectrum for each single theta cycle and normalized to theta phase.

In this case, you need extract instanous theta phase based on hilbert transform. So first filtered the LFP in theta band, and use hilbert transform. This part is not invovled in my uploaded code. An simple implement is following, one can use other filter for sure.

thetaband=[6 12]; %%%%6-12 Hz band passing

samplingRate=1000; %%%%1000Hz;

[b,a]=butter(4,LFP,thetaband/(samplingRate/2),'bandpass');

LFPtheta=filterfilter(b,a,LFPtheta)

thetaPhase=angle(hilbert(LFPtheta));

Then in my paper, each trough of LFP theta was considered as start/end of one single theta cycle. And each single theta cycle was divided into 20 equal phase bins, cross-spectrum was averaged with each phase bin. Then we would have 81(frequcies)*20(phases) matrix the WCS of one theta cycle. Calculate those for all single theta cycles. In the paper, B_k is the cross spectrum, thus W_k=angle(B_k) is the phase. Function angle.m from Mathworks is for that.

Step3: Calculate the PPC.

Calculate PPC for each given phase bin and each given frequeicies based on formula listed in Pair-wise phase consistency for LFP-LFP phase synchrony session in the paper. W_kW_j in the formula is the dot product, where W_kW_j= cos(W_k)*cos(W_j)+sin(W_k)*sin(W_j).

Noted that, the function LU_WCSThetaExtract.txt uploaded combined step1 and step2 above, input variable LFP1, LFP2, LFPtheta and thetaphase is the same length of vector (or time series). Here assuming the sampling time is form [1:length(LFP1)]/samprate.

In the 140 line, a modified mathwork function was used called coher_newVersion2.m, which I didn't upload for copyright reason (I modified the wcoherence.m function from Mathworks into an old fashion to use). One can simply use wcoher.m (old version from mathworks) or wcoherence.m (new version from mathwork). The purpose of the function here is just to calculate cross-spectrum between two siginals, you can also use other wavelet toolbox instead. Noted that, WCS is the cross spectrum, not the wavelet coherence. I recommend to use the old version wcoher.m if you want use my code LU_WCSThetaExtract.m, since it specified the frequency-scale, change the coher_newVersion2 directly to wcoher in line 140 should work.