function [F_mat,bs_begin_sample_final,F,T,SW_definition_thresh_FR,instantaneous_freq_per_bs,no_spike_segments,spikes_per_larger_bs,F_mat_seg,F_seg]=Calculate_Spike_Power_Spectra(spikes_per_ms,ms_begin_sample,binsize,settings,profile_type,bandpass_SW,bandstop_SW,amp_percent_thresh,NREM_bins,segmented_PS);
divide_to_1_s_bins=0;
use_customed_pwelch=0;
end_freq=40; %121518
% This function computes the power spectrum of a spike train. It basically receives the spike count per each ms,
% and the sample corresponding for each ms in the data. It does the parcing to binsizes itself, in a way that
% influences the highest frequency etc. These are specified in binsize and will be explained hereafter.
% "binsize", as opposed to binsize_for_FR_calculation below, is used to set
% the frequency resolution of the power spectra. This is the number of
% seconds I look at to compute the PS. for number of seconds A with
% sampling rate r (set by binsize_for_FR_calculation), I get A*r freq
% points. Entire spectrum is between -(1/r)/2 to (1/r)/2 Hz in the frequency
% domain, i.e. a span of r Hz. A*r freq points in this span give 1/A Hz
% resolution.
if profile_type==1 % sleep. No event locking.
binsize_for_FR_calculation=0.005; % used to be 5. see what happens. in seconds. This binsize determines the amount of 1-ms bins that will be chunked together. Of course, this determines directly the new "sampling rate" of the firing rate vector to be later analyzed. Therefore, here you set the maximal frequency you will see in your spectrum. 100 ms = highest frequency of (1000/100)/2=5 Hz.
if mod(length(spikes_per_ms),1000*binsize_for_FR_calculation)~=0 % first append zeros to create an integer number of bins.
spikes_per_ms_appended=[spikes_per_ms,zeros(1,1000*binsize_for_FR_calculation-mod(length(spikes_per_ms),1000*binsize_for_FR_calculation))]; %complement to nearest binsize*1000 ms multiplier
else
spikes_per_ms_appended=spikes_per_ms;
end
if binsize_for_FR_calculation==0.001
spikes_per_bs=reshape(spikes_per_ms_appended,binsize_for_FR_calculation*1000,length(spikes_per_ms_appended)/(1000*binsize_for_FR_calculation));
else
spikes_per_bs=sum(reshape(spikes_per_ms_appended,binsize_for_FR_calculation*1000,length(spikes_per_ms_appended)/(1000*binsize_for_FR_calculation)));
end
% spikes_per_bs=sum(reshape(spikes_per_ms_appended,binsize_for_FR_calculation*1000,length(spikes_per_ms_appended)/(1000*binsize_for_FR_calculation))); % 060917, move back to non-comment if there is problems, insteas of entire if statement above
if mod(length(spikes_per_ms),1000*binsize_for_FR_calculation)~=0
spikes_per_bs(1,length(spikes_per_bs))=spikes_per_bs(1,length(spikes_per_bs))/(mod(length(spikes_per_ms),(binsize_for_FR_calculation*1000))/(binsize*1000)); % because zeros are appended, the last entry, meaning the last second in this vector doesn't have the actual number of spikes per (binsize)*second, but the number of spikes per x<(binsize*1000) ms. I divide this number by x/(1000*binsize) to have the approximation of the sp/bs FR.
end
spikes_per_bs=spikes_per_bs/binsize_for_FR_calculation; % to move to sp/s
% made comment in 060917, line below
% bs_begin_sample=ms_begin_sample(mod((1:length(ms_begin_sample)),1000*binsize_for_FR_calculation)==1); % all the entries in milli-second_begin_sample that start a full binsize*second rather than 1 ms. also, notice that the first sample in second_begin_sample, just as milli_second_begin_sample, is always StableStartSample
% if statement below, 060917
if binsize_for_FR_calculation==0.001
bs_begin_sample=ms_begin_sample;
else
bs_begin_sample=ms_begin_sample(mod((1:length(ms_begin_sample)),1000*binsize_for_FR_calculation)==1); % all the entries in milli-second_begin_sample that start a full binsize*second rather than 1 ms. also, notice that the first sample in second_begin_sample, just as milli_second_begin_sample, is always StableStartSample
end
% spikes per bs is a vector of spikes by binsize_for_FR_calculation ms bins. bs begin sample is their starts
if settings(1,2)==1
% basically, the next code duplicates parts of the vector to create
% a bigger vector, made of overlapping segments, when the number of
% overlapping samples is set by the user.
binsize_overlap=round(settings(2,2)/100*binsize/binsize_for_FR_calculation); % num of binsize_for_FR_calculation-ms binsizes to overlap
reach_end=0; counter=0;
while reach_end~=1
counter=counter+1;
bs_begin_sample=[bs_begin_sample(1:(counter*binsize/binsize_for_FR_calculation)),bs_begin_sample(counter*binsize/binsize_for_FR_calculation-binsize_overlap+1:counter*binsize/binsize_for_FR_calculation),bs_begin_sample((counter*binsize/binsize_for_FR_calculation+1):end)];
spikes_per_bs=[spikes_per_bs(1:(counter*binsize/binsize_for_FR_calculation)),spikes_per_bs(counter*binsize/binsize_for_FR_calculation-binsize_overlap+1:counter*binsize/binsize_for_FR_calculation),spikes_per_bs((counter*binsize/binsize_for_FR_calculation+1):end)];
if (length(bs_begin_sample)-counter*binsize/binsize_for_FR_calculation)==0
reach_end=1;
add_samples=0;
elseif ((length(bs_begin_sample)-counter*binsize/binsize_for_FR_calculation)<(binsize/binsize_for_FR_calculation))
reach_end=1;
counter=counter+1;
add_samples=counter*binsize/binsize_for_FR_calculation-length(bs_begin_sample);
end
end
spikes_per_bs_appended=[spikes_per_bs,zeros(1,add_samples)];
spikes_per_larger_bs=reshape(spikes_per_bs_appended,binsize/binsize_for_FR_calculation,counter);
bs_begin_sample_appended=[bs_begin_sample,zeros(1,add_samples)];
bs_begin_sample_larger=reshape(bs_begin_sample_appended,binsize/binsize_for_FR_calculation,counter);
bs_begin_sample_final=bs_begin_sample_larger(1,:);
else % next block creates the binsize over number of segment matrix for cases where overlapping segments are not needed. Two cases: first is for when the existing number of bins does not divide by the requested binsize. In this case additional zero bins are complemented, both for the actual spike count data and for the binsize begin samples. Of course, no new "binsizes" are created so no fake bin begin samples are added. The other case is the simplest, where there is no need to append anything.
if mod(length(spikes_per_bs),binsize/binsize_for_FR_calculation)~=0
spikes_per_bs_appended=[spikes_per_bs,zeros(1,binsize/binsize_for_FR_calculation-mod(length(spikes_per_bs),binsize/binsize_for_FR_calculation))]; %complement to nearest binsize*1000 ms multiplier
spikes_per_larger_bs=reshape(spikes_per_bs_appended,binsize/binsize_for_FR_calculation,length(spikes_per_bs_appended)/(binsize/binsize_for_FR_calculation));
bs_begin_sample_appended=[bs_begin_sample,zeros(1,binsize/binsize_for_FR_calculation-mod(length(bs_begin_sample),binsize/binsize_for_FR_calculation))]; %complement to nearest binsize*1000 ms multiplier
bs_begin_sample_larger=reshape(bs_begin_sample_appended,binsize/binsize_for_FR_calculation,length(bs_begin_sample_appended)/(binsize/binsize_for_FR_calculation));
bs_begin_sample_final=bs_begin_sample_larger(1,:);
else
spikes_per_bs_appended=spikes_per_bs;
spikes_per_larger_bs=reshape(spikes_per_bs_appended,binsize/binsize_for_FR_calculation,length(spikes_per_bs_appended)/(binsize/binsize_for_FR_calculation));
bs_begin_sample_appended=bs_begin_sample;
bs_begin_sample_larger=reshape(bs_begin_sample_appended,binsize/binsize_for_FR_calculation,length(bs_begin_sample_appended)/(binsize/binsize_for_FR_calculation));
bs_begin_sample_final=bs_begin_sample_larger(1,:);
end
end
% calcuate the amplitude threshold for SW definition change back!!!
% squashed_only_NREM_spikes_per_larger_bs=reshape(spikes_per_larger_bs(:,NREM_bins),length(NREM_bins)*size(spikes_per_larger_bs,1),1); % changed 061617, used to be all the neuron activity. now only NREM
% spikes_per_bs_appended_only_NREM_for_SW_amp=squashed_only_NREM_spikes_per_larger_bs-nanmean(squashed_only_NREM_spikes_per_larger_bs);
% filtered_NREM_FR_signal_all_unit=BWfilter_designer(double(spikes_per_bs_appended_only_NREM_for_SW_amp), 1/binsize_for_FR_calculation, bandpass_SW, bandstop_SW, 'BP'); % previously [0.2,2], [0.1,2.5]; 35 and 40 for analysis of entire spectrum
% h_filtered__NREM_FR_signal_all_unit=hilbert(filtered_NREM_FR_signal_all_unit); ins_amp_NREM_FR_thresh=abs(h_filtered__NREM_FR_signal_all_unit);
% SW_definition_thresh_FR=prctile(sort(ins_amp_NREM_FR_thresh),amp_percent_thresh);
% before 061617
spikes_per_bs_appended_for_SW_amp=spikes_per_bs_appended-nanmean(spikes_per_bs_appended);
filtered_FR_signal_all_unit=BWfilter_designer(double(spikes_per_bs_appended_for_SW_amp), 1/binsize_for_FR_calculation, bandpass_SW, bandstop_SW, 'BP'); % previously [0.2,2], [0.1,2.5]; 35 and 40 for analysis of entire spectrum
h_filtered_FR_signal_all_unit=hilbert(filtered_FR_signal_all_unit); ins_amp_FR_thresh=abs(h_filtered_FR_signal_all_unit);
SW_definition_thresh_FR=prctile(sort(ins_amp_FR_thresh),amp_percent_thresh);
% so far, only the spike count was calculated. PS calculation begins now.
F_mat=[];
F_mat_seg=[];
% instantaneous_freq_per_bs=nan(10,size(spikes_per_larger_bs,2)); % calculate the instantaneous frequency in each second, for each bin. This is for the supplementary - to show that the instantaneous frequency is not constant.
instantaneous_freq_per_bs=nan(1,size(spikes_per_larger_bs,2)); % before I used to do the calculation per one second such that I had 10 mean. I think it is rather unnecessary.
no_spike_segments=find(sum(spikes_per_larger_bs,1)==0);
for ww=1:size(spikes_per_larger_bs,2)
% eliminate DC
spikes_per_larger_bs(:,ww)=spikes_per_larger_bs(:,ww)-nanmean(spikes_per_larger_bs(:,ww)); % subtract mean to avoid DC
% for when binsize=2, I want to make it 3 by ZP. It is not
% recommended to use binsize 1 or less. The smaller binsize I think
% is manageable, 2 sec, I zero pad with extra zeros to make 3-s
% segments to get 0.33 Hz resolution, which is minimal. If the
% initial segment is bigger, I do ZP only for the closest power of
% 2 for ease of computation.
if binsize==2
[PXX,F]=periodogram(spikes_per_larger_bs(:,ww),hamming(length(spikes_per_larger_bs(:,ww))),2^(nextpow2(1.5*length(spikes_per_larger_bs(:,ww)))),1/(binsize_for_FR_calculation));
else
if divide_to_1_s_bins
PXX=[];
for b=1:10
a=spikes_per_larger_bs((b-1)*200+1:b*200,ww);
a=a-nanmean(a);
[px,F]=periodogram(a,hamming(length(a)),2^(nextpow2(length(a))),200);
PXX=[PXX,px];
end
PXX=PXX(1:find(F>end_freq,1,'first'),:);
PXX=PXX./repmat(sum(PXX,1),size(PXX,1),1);
PXX=nanmean(PXX,2);
elseif use_customed_pwelch
[PXX,F]=pwelch(spikes_per_larger_bs(:,ww),hamming(length(spikes_per_larger_bs(:,ww))),0.5*length(spikes_per_larger_bs(:,ww)),[0:0.1:end_freq],200);
PXX=PXX';
else
[PXX,F]=periodogram(spikes_per_larger_bs(:,ww),hamming(length(spikes_per_larger_bs(:,ww))),2^(nextpow2(length(spikes_per_larger_bs(:,ww)))),1/(binsize_for_FR_calculation));
if segmented_PS
ZZ=spikes_per_larger_bs(:,ww)';
seglen=length(ZZ)/5;
for kkk=1:5
ZZZ=ZZ(1,(kkk-1)*seglen+1:seglen*kkk);
ZZZ=[ZZZ,zeros(1,4*seglen)];
[PXX_SEG{1,kkk},F_seg]=periodogram(ZZZ-nanmean(ZZZ),hamming(length(ZZZ)),2^(nextpow2(length(ZZZ))),1/(binsize_for_FR_calculation));
PXX_SEG{1,kkk}=PXX_SEG{1,kkk}(1:find(F_seg>end_freq,1,'first'),:);
end
end
end
% [PXX,F]=pwelch(spikes_per_larger_bs(:,ww),hamming(length(spikes_per_larger_bs(:,ww))),0,2^(nextpow2(length(spikes_per_larger_bs(:,ww)))),1/(binsize_for_FR_calculation));
% seg_length=length(spikes_per_larger_bs(:,ww))/2;
% [PXX,F]=pwelch(spikes_per_larger_bs(:,ww),blackman(seg_length),seg_length/5,2^(nextpow2(seg_length)),1/(binsize_for_FR_calculation));
end
% normalize power spectra to get relative power
PXX=PXX/sum(PXX);
F_mat=[F_mat,PXX];
if segmented_PS
for kkk=1:5
PXX_SEG{1,kkk}=PXX_SEG{1,kkk}/sum(PXX_SEG{1,kkk});
end
F_mat_seg=cat(3,F_mat_seg,[PXX_SEG{1,1},PXX_SEG{1,2},PXX_SEG{1,3},PXX_SEG{1,4},PXX_SEG{1,5}]);
end
filtered=BWfilter_designer(double(spikes_per_larger_bs(:,ww)),1/binsize_for_FR_calculation, bandpass_SW, bandstop_SW, 'BP');
h_filtered=hilbert(filtered');
h_amp=abs(h_filtered);
ins_freq_h_filtered=[0,(1/binsize_for_FR_calculation)/(2*pi)*diff(unwrap(angle(h_filtered)))];
ins_freq_h_filtered(find(ins_freq_h_filtered<=0 | ins_freq_h_filtered>=bandpass_SW(1,2)))=NaN; %changed from 1 to bandpass_SW(1,2) 121418
% instantaneous_freq_per_bs(:,ww)=nanmean(reshape(ins_freq_h_filtered,length(ins_freq_h_filtered)/10,10),1)'; % This has to do with the 10-mean thing above. not necessary.
instantaneous_freq_per_bs(:,ww)=nanmean(ins_freq_h_filtered(find(h_amp>SW_definition_thresh_FR))); % take only the instantaneous frequency for samples that are bigger than the thresh defined above for the entire oscillation
% calculate instantaneous frequency
end
% F_mat=10*log10(abs(real(F_mat))); % convert to dB. No need to convert
% to dB as I work with relative power - every freq is given a percent
% representing its part of the total power.
T=0;
if ~segmented_PS
F_mat_seg=[]; F_seg=[];
end
elseif profile_type==2 % Event. Event locking.
spikes_per_ms=1000*reshape(spikes_per_ms',1,size(spikes_per_ms,1)*size(spikes_per_ms,2));
smoothed_spikes_per_ms=filtfilt(gausswin(settings.kernel_length)/sum(gausswin(settings.kernel_length)),1,spikes_per_ms);
smoothed_spikes_per_ms=smoothed_spikes_per_ms-nanmean(smoothed_spikes_per_ms);
[S,F,T,F_mat]=spectrogram(smoothed_spikes_per_ms, hann(2^nextpow2(binsize)), 2^nextpow2(binsize)/2, 2^nextpow2(binsize), 1000, 'yaxis');
T=T-settings.PRESEC;
bs_begin_sample_final=0;
end
end