s

Lindsey_WF_phasefinder.m

@@ -0,0 +1,114 @@
+%inputLFPep input peak theta freq      
+Theta_filt = designfilt('bandpassiir','FilterOrder',12, ...
+    'HalfPowerFrequency1',thetafreq-2, 'HalfPowerFrequency2',thetafreq+2, ...
+    'SampleRate',sFreq,'designmethod', 'butter');
+
+narrow_filt = designfilt('bandpassiir','FilterOrder',12, ...
+    'HalfPowerFrequency1',thetafreq-(2*1.2), 'HalfPowerFrequency2',thetafreq+(2*1.2), ... making 20% bigger
+    'SampleRate',sFreq,'designmethod', 'butter');
+
+broad_filt = designfilt('bandpassiir','FilterOrder',12, ...
+    'HalfPowerFrequency1',1, 'HalfPowerFrequency2',28, ...%based on bellusco
+    'SampleRate',sFreq,'designmethod', 'butter');
+
+narrow(:,1)=LFPep(:,1);
+broad(:,1)=LFPep(:,1);
+    narrow(:,2)=filtfilt(narrow_filt,LFPep(:,2));
+        broad(:,2)=filtfilt(broad_filt,LFPep(:,2));
+
+
+    theta_phase=[];
+    theta_phase(:,1)=LFPep(:,1);
+
+    theta_phase(:,2)=angle(hilbert(filtfilt(Theta_filt,LFPep(:,2))));
+
+%PEAK BASED ON FITLERED
+[PeaksIdx,TroughsIdx,~, ~] = Find_peaks_troughs_zeros(narrow(:,2));
+[~,~,zupix, zdownix] = Find_peaks_troughs_zeros(broad(:,2));
+
+figure;plot(LFPep(:,1),LFPep(:,2))
+testix=LFPep(:,1)>3.97e9 & LFPep(:,1)<3.9735e9;
+
+figure;plot(LFPep(testix,1),LFPep(testix,2));
+testLFP=LFPep(testix,:);
+testnarrow=narrow(testix,:);
+testbroad=broad(testix,:);
+
+[PeaksIdx,TroughsIdx,~, ~] = Find_peaks_troughs_zeros(testnarrow(:,2));
+ [pkix,pks]= findpeaks(testnarrow(:,2),testnarrow(:,1)); %now locs is expressed in time
+ [pkvals,pkix]= findpeaks(testnarrow(:,2)); %now locs is expressed in time
+
+
+ [trix,trg]= findpeaks(-(testnarrow(:,2)),testnarrow(:,1)); %now locs is expressed in time
+ [trvals,trix]= findpeaks(-(testnarrow(:,2))); %now locs is expressed in time
+
+%narrow now just needs to back the peaks and troughs up to the nearest
+%actual peak in broad
+
+[~,~,zupix, zdownix] = Find_peaks_troughs_zeros(testbroad(:,2));
+figure
+plot(testLFP(:,1),testLFP(:,2))
+hold on
+plot(testnarrow(:,1),testnarrow(:,2),'k')
+plot(testbroad(:,1),testbroad(:,2),'g')
+
+plot(testLFP(pkix,1),testLFP(pkix,2),'ob')
+plot(testLFP(pkix,1),testLFP(trix,2),'ob')
+plot(testLFP(zupix,1),testLFP(zupix,2),'or')
+plot(testLFP(zdownix,1),testLFP(zdownix,2),'or')
+
+                [~,locs]= findpeaks(narrow(:,2),narrow(:,1)); %now locs is expressed in time
+                centerpeak=ceil(length(locs)/2);
+                theta_peak=locs(centerpeak);
+                theta_peakm1=locs(centerpeak-1);
+                theta_peakp1=locs(centerpeak+1);
+
+        %trough  1 and 2 based on theta   
+                [~,locs2]= findpeaks(-(filtered_theta(time_ix,2)),LFPep(time_ix,1)); %now locs is expressed in time
+                backix=locs2<theta_peak;
+                theta_trough=max(locs2(backix));
+                frontix=locs2>theta_peak;
+                theta_trough2=min(locs2(frontix));
+
+       %ninety degress
+            snip= Restrict(broad_filtered,theta_trough,theta_peak);
+             ninety_ix=dsearchn(snip(:,2),0);
+            ninety_time=snip(ninety_ix,1);
+
+       %180 degrees
+       %find second trough, then same as for 90
+                snip2= Restrict(broad_filtered,theta_peak,theta_trough2);
+             oneeighty_ix=dsearchn(snip2(:,2),0);
+            oneeighty_time=snip2(oneeighty_ix,1);
+
+       %zero degrees
+       snip3=Restrict(broad_filtered,theta_peakm1,theta_trough);
+              zero_ix=dsearchn(snip3(:,2),0);
+            zero_time=snip3(zero_ix,1);
+
+         %two70 degrees
+       snip4=Restrict(broad_filtered,theta_trough2,theta_peakp1);
+              two70_ix=dsearchn(snip4(:,2),0);
+            two70_time=snip4(two70_ix,1);          
+
+
+            % NOW REAL PEAK BASED ON MAX and MIN between zero crossings
+            Wave_area=Restrict(broad_filtered,ninety_time,oneeighty_time);
+            [val,ix]=max(Wave_area(:,2));
+            time_peak=Wave_area(ix,1);
+
+            %real trough 1 base on zero and 90
+            Wave_area=Restrict(broad_filtered,zero_time,ninety_time);
+            [val,ix]=min(Wave_area(:,2));
+            time_trough=Wave_area(ix,1);
+
+            %real trough 2
+            Wave_area=Restrict(broad_filtered,oneeighty_time,two70_time);
+            [val,ix]=min(Wave_area(:,2));
+            time_trough2=Wave_area(ix,1);    
+
+
+             Quad1_us= ninety_time-time_trough  ; %Quad 1: Trough 1 to ninety
+            Qaud2_us= time_peak-ninety_time ;   %Quad 2: Ninety to peak
+            Quad3_us= oneeighty_time-time_peak ;  %Quad 3: Peak to 180
+             Quad4_us=time_trough2-oneeighty_time;   %Quad 4: 180 to trough 2

Phase_pow_fq_detector_waveshape.m

@@ -0,0 +1,175 @@
+function [phase, pk_ix, tr_ix, assym, updown_ix] = ...
+    Phase_pow_fq_detector_waveshape(wideLFP,narrowLFP,sFreq)
+% function [phase, pk_ix, tr_ix, assym, updown_ix] = ...
+%     Phase_pow_fq_detector_waveshape(wideLFP,narrowLFP,sFreq)
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% Make narrowLFP not too narrow. perhaps 20% wider than your typical range
+% for the frequency. WideLFP - get rid of slow changes in mean - perhaps
+% medfilt before sending it in.
+%
+% NOTE: Upsample this if you want better time/phase resolution.
+%
+% Make wide not complete wide band. Perhaps 200% wider than the target
+% frequency.
+%
+% SEE BULLSCIO BUZSAKI PAPER
+% Belluscio, M. A., Mizuseki, K., Schmidt, R., Kempter, R., & Buzsáki, G. (2012). Cross-frequency phase-phase coupling between ? and ? oscillations in the hippocampus. The Journal of Neuroscience : The Official Journal of the Society for Neuroscience, 32(2), 423–435. https://doi.org/10.1523/JNEUROSCI.4122-11.2012
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% Cowen 2020 - updated to uncorporate zero crossings.
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+if nargin ==0
+    % Validate with some real data...
+    tst = load('I:\waveformstuff.mat');
+    figure
+    plot(tst.LFPep(:,1),tst.LFPep(:,2))
+    hold on
+    plot(tst.narrow(:,1),tst.narrow(:,2))
+    plot(tst.broad(:,1),tst.broad(:,2))
+    sFreq = 1e6/median(diff(tst.broad(:,1)));
+    wideLFP = tst.broad(:,2);
+    narrowLFP = tst.narrow(:,2);
+end
+
+assym = []; inst_fq2 = [];
+if nargin < 3
+    sFreq = 1;
+end
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+[~,~,zup_ix,zdown_ix] = Find_peaks_troughs_zeros(narrowLFP);
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+new_zdown_ix = zeros(size(zup_ix));
+for ii = 1:length(zup_ix)-1
+    tmp = find(zdown_ix > zup_ix(ii),1,'first');
+    if isempty(tmp)
+        break
+    else
+        new_zdown_ix(ii) = zdown_ix(tmp);
+    end
+end
+zup_ix = zup_ix(1:(ii-1));
+zdown_ix = new_zdown_ix(1:(ii-1));
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% always start with the peak - the rise up, end with a trough. Keeps things consistent.
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+if zup_ix(1) > zdown_ix(1)
+    zdown_ix(1) = [];
+end
+if zup_ix(end) > zdown_ix(end)
+    zup_ix(end) = [];
+end
+updown_ix = sort([zup_ix;zdown_ix]);
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% Now deal with the WIDE filtered data. Find peaks and troughs...
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% Find the zero crossings for the wide data...
+[~,~,zup_ix_wide,zdown_ix_wide] = Find_peaks_troughs_zeros(wideLFP);
+% find the closes zup_ix_wide for each zup_ix...
+tmp = ones(length(zup_ix),1);
+for ii = 1:length(zup_ix)
+    ix = binsearch(zup_ix_wide,zup_ix(ii));
+    tmp(ii) = zup_ix_wide(ix);
+end
+zup_wide_ix = tmp;
+
+tmp = ones(length(zdown_ix),1);
+for ii = 1:length(zdown_ix)
+    ix = binsearch(zdown_ix_wide,zdown_ix(ii));
+    tmp(ii) = zdown_ix_wide(ix);
+end
+zdown_ix_wide = tmp;
+
+updown_wide_ix = sort([zup_wide_ix;zdown_ix_wide]);
+
+cnt = 1;
+pk_ix = nan(round(length(updown_ix)/2),1);
+tr_ix = nan(round(length(updown_ix)/2),1);
+for ii = 1:2:length(updown_ix)-2
+    [~, m] = max(wideLFP(updown_ix(ii):updown_ix(ii+1)));
+    pk_ix(cnt) = m + updown_ix(ii) - 1;
+    [~, m] = min(wideLFP(updown_ix(ii+1):updown_ix(ii+2)));
+    tr_ix(cnt) = m + updown_ix(ii+1) - 1;
+    cnt = cnt + 1;
+end
+tr_ix = tr_ix(1:(cnt-1));
+pk_ix = pk_ix(1:(cnt-1));
+fall = (tr_ix - pk_ix);
+rise = (pk_ix(2:end) - tr_ix(1:end-1));
+rise = [rise(1); rise];
+% get rid of strange situations where the the peak and trough are the same
+% value. 
+BIX = fall <= 0 | rise <=0;
+pk_ix = pk_ix(~BIX);
+tr_ix = tr_ix(~BIX);
+fall = (tr_ix - pk_ix);
+rise = (pk_ix(2:end) - tr_ix(1:end-1));
+rise = [rise(1); rise];
+
+assym = [pk_ix (fall)./(rise + fall) - .5] ;
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%% Problems - at every 90 degrees.
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+phase_old = nan(size(wideLFP));
+phase = nan(size(wideLFP));
+for ii = 1:length(pk_ix)
+    % find the zero crossing after the peak.
+    zcix = find(zdown_ix_wide>pk_ix(ii),1,'first');
+    zc = zdown_ix_wide(zcix);
+    if zc == tr_ix(ii)
+        zc = tr_ix(ii) - round((tr_ix(ii) - pk_ix(ii))/2);
+    end
+
+    ix = pk_ix(ii):tr_ix(ii);
+    phase(ix) = interp1([pk_ix(ii) zc tr_ix(ii)],[0 pi/2 pi],ix(:),'linear');
+
+    phase_old(ix) = linspace(0,pi,length(ix));
+end
+for ii = 1:length(tr_ix)-1
+    zcix = find(zup_ix_wide>tr_ix(ii),1,'first');
+    zc = zup_ix_wide(zcix);
+    if zc == pk_ix(ii+1)
+        zc = pk_ix(ii+1) - round((pk_ix(ii+1) - tr_ix(ii))/2);
+    end
+    w
+    ix = tr_ix(ii):pk_ix(ii+1);
+
+    phase(ix) = interp1([tr_ix(ii) zc pk_ix(ii+1)],[pi 1.5*pi 2*pi],ix,'linear');
+
+    phase_old(ix) = linspace(pi,2*pi,length(ix));
+end
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+if 0
+    figure 
+    % NOTE: THe clustering around 0, 180 degrees makes sense as we are
+    % using linspace to assign phase which GUARANTEES that you will have a
+    % point at 0 180 for EVERY cycle. The other phases will be less
+    % represented because linspace will assign different phases depending
+    % on the number of points. Could up the sampling rate and this would
+    % reduce the problem a little.
+    histogram(phase,700)
+
+    figure
+    histogram(assym(:,2),70)
+
+    %
+    figure
+    plot(wideLFP,'k')
+    hold on
+    plot(wideLFP,'k.')
+
+    plot(narrowLFP,'b')
+    plot(updown_ix(1:2:end),narrowLFP(updown_ix(1:2:end)),'ro')
+    plot(updown_ix(2:2:end),narrowLFP(updown_ix(2:2:end)),'b*')
+    plot(updown_wide_ix(1:2:end),wideLFP(updown_wide_ix(1:2:end)),'mp')
+    plot(updown_wide_ix(2:2:end),wideLFP(updown_wide_ix(2:2:end)),'m+')
+    plot_horiz_line_at_zero()
+
+    plot(tr_ix,wideLFP(tr_ix),'m^')
+    plot(pk_ix,wideLFP(pk_ix),'g^')
+    set(gca,'Ylim',prctile(wideLFP(1:40:end),[.1 99.9]))
+
+    plot(rad2deg(phase),'r-')
+    hold on
+    plot(rad2deg(phase_old),'b-')
+    %
+end

Prep_EEG_4_PAC.m

@@ -0,0 +1,26 @@
+function [pwr,phase]=Prep_EEG_4_PAC(sFreq,freq4power,freq4phase,EEG)
+
+% sFreq = 1000;
+time = -1:1/sFreq:1;
+half_of_wavelet_size = (length(time)-1)/2;
+n_wavelet     = length(time);
+n_data        = length(EEG);
+n_convolution = n_wavelet+n_data-1;
+fft_data = fft(EEG,n_convolution);
+% freq4phase=8;
+% freq4power=30;
+% wavelet for phase and its FFT
+wavelet4phase = exp(2*1i*pi*freq4phase.*time) .* exp(-time.^2./(2*(4/(2*pi*freq4phase))^2));
+fft_wavelet4phase = fft(wavelet4phase,n_convolution);
+
+% wavelet for power and its FFT
+wavelet4power = exp(2*1i*pi*freq4power.*time) .* exp(-time.^2./(2*(4/(2*pi*freq4power))^2));
+fft_wavelet4power = fft(wavelet4power,n_convolution);
+
+% get phase values
+convolution_result_fft = ifft(fft_wavelet4phase.*fft_data,n_convolution);
+phase = angle(convolution_result_fft(half_of_wavelet_size+1:end-half_of_wavelet_size));
+
+% get power values (note: 'power' is a built-in function so we'll name this variable 'amp')
+convolution_result_fft = ifft(fft_wavelet4power.*fft_data,n_convolution);
+pwr = abs(convolution_result_fft(half_of_wavelet_size+1:end-half_of_wavelet_size)).^2;