CALCPOWER - Compute spectral analysis of DynaSim data
Usage:
data = dsCalcPower(data,'option',value)
Inputs:
- data: DynaSim data structure (see dsCheckData)
- options:
'variable' : name of field containing data on which to
calculate firing rates (default: *_spikes or
first variable in data.labels)
'time_limits' : [beg,end] (units of data.time)
'smooth_factor' : number of samples for smoothing the spectrum (default: 5)
- tip: set to 1 to avoid smoothing
- options for peak detection:
'min_peak_frequency' : Hz, min frequency for peak detection (default: 2)
'max_peak_frequency' : Hz, max frequency for peak detection (default: 150)
'peak_threshold_prctile': percentile for setting power threshold for peak
detection (default: 95)
'peak_area_width' : Hz, size of frequency bin (centered on peak)
over which to calculate area under spectrum (default: 5)
'exclude_data_flag' : whether to remove simulated data from result
structure (default: 0)
Outputs:
- data structure organization:
data.VARIABLE_Power_SUA.frequency: TODO
data.VARIABLE_Power_SUA.PeakArea: area under spectrum around peak (one value per cell)
data.VARIABLE_Power_SUA.PeakFreq: frequency of spectral power (one value per cell)
data.VARIABLE_Power_SUA.Pxx: spectral power
- if populations present, data also includes:
data.VARIABLE_Power_MUA.frequency: TODO
data.VARIABLE_Power_MUA.PeakArea: TODO
data.VARIABLE_Power_MUA.PeakFreq: TODO
data.VARIABLE_Power_MUA.Pxx: spectrum of the mean waveform
- population mean spectrum of the individual waveforms can be
calculated using "mean(data.VARIABLE_Power_MUA.Pxx,2)".
- Note:
- "VARIABLE" can be specified as the name of a variable listed in
data.labels, a cell array of string listing variable names, or as a
regular expression pattern for identifying variables to process. See
dsSelectVariables for more info on supported specifications.
Examples:
s=[];
s.populations(1).name='E';
s.populations(1).equations='dv[2]/dt=@current+10; {iNa,iK}; v(0)=-65';
s.populations(2).name='I';
s.populations(2).equations='dv/dt=@current+10; {iNa,iK}; v(0)=-65';
data=dsSimulate(s,'tspan',[0 1000]);
data=dsCalcPower(data,'variable','v');
% Plot the spectrum of the E-cell average population voltage
figure; plot(data.E_v_Power_MUA.frequency,data.E_v_Power_MUA.Pxx);
xlabel('frequency (Hz)'); ylabel('power'); xlim([0 200]);
See also: PlotPower, dsAnalyzeStudy, dsSimulate, dsCheckData, dsSelectVariables
Author: Jason Sherfey, PhD <jssherfey@gmail.com>
Copyright (C) 2016 Jason Sherfey, Boston University, USA0001 function data = dsCalcPower(data, varargin) 0002 %CALCPOWER - Compute spectral analysis of DynaSim data 0003 % 0004 % Usage: 0005 % data = dsCalcPower(data,'option',value) 0006 % 0007 % Inputs: 0008 % - data: DynaSim data structure (see dsCheckData) 0009 % - options: 0010 % 'variable' : name of field containing data on which to 0011 % calculate firing rates (default: *_spikes or 0012 % first variable in data.labels) 0013 % 'time_limits' : [beg,end] (units of data.time) 0014 % 'smooth_factor' : number of samples for smoothing the spectrum (default: 5) 0015 % - tip: set to 1 to avoid smoothing 0016 % - options for peak detection: 0017 % 'min_peak_frequency' : Hz, min frequency for peak detection (default: 2) 0018 % 'max_peak_frequency' : Hz, max frequency for peak detection (default: 150) 0019 % 'peak_threshold_prctile': percentile for setting power threshold for peak 0020 % detection (default: 95) 0021 % 'peak_area_width' : Hz, size of frequency bin (centered on peak) 0022 % over which to calculate area under spectrum (default: 5) 0023 % 'exclude_data_flag' : whether to remove simulated data from result 0024 % structure (default: 0) 0025 % 0026 % Outputs: 0027 % - data structure organization: 0028 % data.VARIABLE_Power_SUA.frequency: TODO 0029 % data.VARIABLE_Power_SUA.PeakArea: area under spectrum around peak (one value per cell) 0030 % data.VARIABLE_Power_SUA.PeakFreq: frequency of spectral power (one value per cell) 0031 % data.VARIABLE_Power_SUA.Pxx: spectral power 0032 % - if populations present, data also includes: 0033 % data.VARIABLE_Power_MUA.frequency: TODO 0034 % data.VARIABLE_Power_MUA.PeakArea: TODO 0035 % data.VARIABLE_Power_MUA.PeakFreq: TODO 0036 % data.VARIABLE_Power_MUA.Pxx: spectrum of the mean waveform 0037 % - population mean spectrum of the individual waveforms can be 0038 % calculated using "mean(data.VARIABLE_Power_MUA.Pxx,2)". 0039 % - Note: 0040 % - "VARIABLE" can be specified as the name of a variable listed in 0041 % data.labels, a cell array of string listing variable names, or as a 0042 % regular expression pattern for identifying variables to process. See 0043 % dsSelectVariables for more info on supported specifications. 0044 % 0045 % Examples: 0046 % s=[]; 0047 % s.populations(1).name='E'; 0048 % s.populations(1).equations='dv[2]/dt=@current+10; {iNa,iK}; v(0)=-65'; 0049 % s.populations(2).name='I'; 0050 % s.populations(2).equations='dv/dt=@current+10; {iNa,iK}; v(0)=-65'; 0051 % data=dsSimulate(s,'tspan',[0 1000]); 0052 % data=dsCalcPower(data,'variable','v'); 0053 % % Plot the spectrum of the E-cell average population voltage 0054 % figure; plot(data.E_v_Power_MUA.frequency,data.E_v_Power_MUA.Pxx); 0055 % xlabel('frequency (Hz)'); ylabel('power'); xlim([0 200]); 0056 % 0057 % See also: PlotPower, dsAnalyzeStudy, dsSimulate, dsCheckData, dsSelectVariables 0058 % 0059 % Author: Jason Sherfey, PhD <jssherfey@gmail.com> 0060 % Copyright (C) 2016 Jason Sherfey, Boston University, USA 0061 0062 %% 1.0 Check inputs 0063 options=dsCheckOptions(varargin,{... 0064 'variable',[],[],... 0065 'time_limits',[-inf inf],[],... 0066 'smooth_factor',5,[],... % number of samples for smoothing the spectrum 0067 'min_peak_frequency',1,[],... % Hz, min frequency for peak detection 0068 'max_peak_frequency',200,[],... % Hz, max frequency for peak detection 0069 'peak_threshold_prctile',95,[],... % percentile for setting power threshold for peak detection 0070 'peak_area_width',5,[],... % Hz, size of frequency bin (centered on peak) over which to calculate area under spectrum 0071 'exclude_data_flag',0,{0,1},... 0072 'timeBandwidthProduct',[],[],... % time-bandwidth product for multi-taper method 0073 'output_suffix','',[],... 0074 'auto_gen_test_data_flag',0,{0,1},... 0075 },false); 0076 0077 %% auto_gen_test_data_flag argin 0078 if options.auto_gen_test_data_flag 0079 varargs = varargin; 0080 varargs{find(strcmp(varargs, 'auto_gen_test_data_flag'))+1} = 0; 0081 varargs(end+1:end+2) = {'unit_test_flag',1}; 0082 argin = [{data}, varargs]; % specific to this function 0083 end 0084 0085 data = dsCheckData(data, varargin{:}); 0086 % note: calling dsCheckData() at beginning enables analysis function to 0087 % accept data matrix [time x cells] in addition to DynaSim data structure. 0088 0089 if numel(data)>1 0090 % use dsAnalyzeStudy to recursively call dsCalcPower on each data set 0091 data=dsAnalyzeStudy(data,@dsCalcPower,varargin{:}); 0092 return; 0093 end 0094 0095 % time parameters 0096 time = data.time; % time vector 0097 dt = time(2)-time(1); % time step 0098 % ntime=length(time); % number of time points in full data set 0099 t1=nearest(time,options.time_limits(1)); % index to first sample 0100 t2=nearest(time,options.time_limits(2)); % index to last sample 0101 nsamp=t2-t1+1; % number of samples for spectral estimate 0102 0103 % frequency parameters 0104 Fs = fix(1/(dt/1000)); % effective sampling frequency 0105 Fmin=options.min_peak_frequency; % min frequency for peak detection 0106 Fmax=options.max_peak_frequency; % max frequency for peak detection 0107 thresh_prctile=options.peak_threshold_prctile; % percentile for setting power threshold for peak detection 0108 smooth_factor=options.smooth_factor; % number of samples to smooth spectrum 0109 Fwin=options.peak_area_width; % size of frequency bin (centered on peak) for calculating area under spectrum 0110 FreqRange=[max(Fmin,2/time(end)) Fmax]; % range to look for spectral peaks 0111 NFFT=2^(nextpow2(nsamp-1)-1);%2); % <-- use higher resolution to capture STO freq variation 0112 % WINDOW=2^(nextpow2(NFFT-1)-3); 0113 % NOVERLAP=[]; % spectral parameters 0114 NW = options.timeBandwidthProduct; 0115 0116 %% 2.0 set list of variables to process as cell array of strings 0117 options.variable=dsSelectVariables(data(1).labels,options.variable, varargin{:}); 0118 0119 %% 3.0 calculate power spectrum for each variable 0120 if ~isfield(data,'results') 0121 data.results={}; 0122 end 0123 0124 warning off 0125 for v=1:length(options.variable) 0126 % extract this data set 0127 var=options.variable{v}; 0128 dat=data.(var); 0129 0130 % determine how many cells are in this data set 0131 ncells=size(dat,2); 0132 0133 % preallocation 0134 PeakFreq=nan(1,ncells); 0135 PeakArea=nan(1,ncells); 0136 0137 % SUA spectra: loop over cells 0138 for i=1:ncells 0139 % select data 0140 X=detrend(dat(t1:t2,i)); % detrend the data 0141 % calculate spectral estimate 0142 if strcmp(reportUI,'matlab') 0143 [tmpPxx,f] = pmtm(X, NW, NFFT, Fs); % calculate power 0144 elseif exist('pwelch') == 2 % 'pwelch is in Octave's path 0145 [tmpPxx,f] = pwelch(X,NFFT,[],NFFT,Fs); % calculate power in octave (pmtm is not implemented yet) 0146 elseif exist('pwelch') ~= 2 % 'pwelch is not in Octave's path 0147 try 0148 pkg load signal; % trying to load octave forge 'signal' package before using pwelch function 0149 fprintf('''pmtm'' function for spectral analysis not available in Octave, using pwelch.\n') 0150 [tmpPxx,f] = pwelch(X,NFFT,[],NFFT,Fs); % calculate power in octave (pmtm is not implemented yet) 0151 catch 0152 error('pwelch function is needed for spectral analysis in Octave, please install the signal package from Octave Forge'); 0153 end 0154 end 0155 0156 if i==1 0157 % get size of spectrum and preallocate result matrix 0158 nfreq=length(f); 0159 Pxx=nan(nfreq,ncells); 0160 end 0161 0162 if all(isnan(tmpPxx(:))) 0163 tmpPxx=zeros(size(tmpPxx)); 0164 end 0165 0166 if ~isa(tmpPxx,'double') 0167 % convert to double precision 0168 tmpPxx=double(tmpPxx); 0169 end 0170 0171 % smooth the spectrum 0172 if smooth_factor>1 && strcmp(reportUI,'matlab') 0173 tmpPxx=smooth(tmpPxx,smooth_factor); 0174 else 0175 tmpPxx=lsmooth(tmpPxx,smooth_factor); 0176 end 0177 0178 % Peak Detection: 0179 % select range of frequencies over which to look for peaks 0180 sel = find(FreqRange(1)<=f & f<=FreqRange(end)); 0181 0182 % set threshold for peak detection 0183 ht=prctile(tmpPxx(sel),thresh_prctile); % ht=prctile(log10(tmpPxx(sel)),thresh_prctile); 0184 0185 if ~isnan(ht) 0186 % get index of peaks in range over threshold 0187 if strcmp(reportUI,'matlab') 0188 [linPeakPower,PPind]=findpeaks(tmpPxx(sel),'MinPeakHeight',ht,'NPeaks',3); % [PeakPower,PPind]=findpeaks(log10(tmpPxx(sel)),'MinPeakHeight',ht,'NPeaks',3); 0189 else 0190 [linPeakPower,PPind]=findpeaks(tmpPxx(sel),'MinPeakHeight',ht,'MinPeakDistance',0,'MinPeakWidth',0); 0191 end 0192 PeakPower = log10(linPeakPower); 0193 else 0194 PPind=[]; 0195 end 0196 0197 if ~isempty(PPind) 0198 % if multiple peaks, only consider the largest 0199 if numel(PPind)>1 0200 PPind=PPind(max(PeakPower)==PeakPower); %PPind=PPind(1); 0201 end 0202 0203 % get frequency at that index 0204 PeakFreq(i) = f(sel(PPind)); 0205 0206 % set limits for calculating area under spectrum around peak 0207 flo=PeakFreq(i)-Fwin/2; 0208 fhi=PeakFreq(i)+Fwin/2; 0209 sel2=(flo<=f & f<=fhi); 0210 % calculate area under spectrum around peak 0211 PeakArea(i) = sum(tmpPxx(sel2))*(f(2)-f(1)); 0212 else 0213 PeakFreq(i)=nan; 0214 PeakArea(i)=nan; 0215 end 0216 % Store results 0217 Pxx(:,i)=tmpPxx; 0218 end 0219 % ----------------------------------------------------- 0220 % Repeat spectral estimate for MUA: 0221 if ncells==1 0222 % same as SUA 0223 Pxx_mean=Pxx; 0224 Pxx_mean_PeakFreq=PeakFreq; 0225 Pxx_mean_PeakArea=PeakArea; 0226 else 0227 if ~strcmp(reportUI,'matlab') && exist('nanmean') ~= 2 % 'nanmean is not in Octave's path 0228 try 0229 pkg load statistics; % trying to load octave forge 'statistics' package before using nanmean function 0230 catch 0231 error('nanmean function is needed, please install the statistics package from Octave Forge'); 0232 end 0233 end 0234 % calculate MUA 0235 X=detrend(nanmean(dat(t1:t2,:),2)); % detrend the data 0236 0237 % calculate spectral estimate 0238 if ~strcmp(reportUI,'matlab') && exist('pwelch') ~= 2 % 'pwelch is not in Octave's path 0239 try 0240 pkg load signal; % trying to load octave forge 'signal' package before using pwelch function 0241 catch 0242 error('pwelch function is needed for spectral analysis in Octave, please install the signal package from Octave Forge'); 0243 end 0244 end 0245 [tmpPxx,f] = pwelch(X,NFFT,[],NFFT,Fs); % calculate power 0246 if all(isnan(tmpPxx(:))) 0247 tmpPxx=zeros(size(tmpPxx)); 0248 end 0249 0250 if ~isa(tmpPxx,'double') 0251 % convert to double precision 0252 tmpPxx=double(tmpPxx); 0253 end 0254 0255 % smooth the spectrum 0256 if smooth_factor>1 && strcmp(reportUI,'matlab') 0257 tmpPxx=smooth(tmpPxx,smooth_factor); 0258 else 0259 tmpPxx=lsmooth(tmpPxx,smooth_factor); 0260 end 0261 0262 % Peak Detection: 0263 % select range of frequencies over which to look for peaks 0264 sel = find(FreqRange(1)<=f & f<=FreqRange(end)); 0265 0266 % set threshold for peak detection 0267 ht=prctile(tmpPxx(sel),thresh_prctile); % ht=prctile(log10(tmpPxx(sel)),thresh_prctile); 0268 if ~isnan(ht) 0269 % get index of peaks in range over threshold 0270 if strcmp(reportUI,'matlab') 0271 [linPeakPower,PPind]=findpeaks(tmpPxx(sel),'MinPeakHeight',ht,'NPeaks',3); % [PeakPower,PPind]=findpeaks(log10(tmpPxx(sel)),'MinPeakHeight',ht,'NPeaks',3); 0272 else 0273 [linPeakPower,PPind]=findpeaks(tmpPxx(sel),'MinPeakHeight',ht,'MinPeakDistance',0,'MinPeakWidth',0); 0274 end 0275 PeakPower = log10(linPeakPower); 0276 else 0277 PPind=[]; 0278 end 0279 0280 if ~isempty(PPind) 0281 % if multiple peaks, only consider the largest 0282 if numel(PPind)>1 0283 PPind=PPind(max(PeakPower)==PeakPower); %PPind=PPind(1); 0284 end 0285 0286 % get frequency at that index 0287 Pxx_mean_PeakFreq = f(sel(PPind)); 0288 0289 % set limits for calculating area under spectrum around peak 0290 flo=Pxx_mean_PeakFreq-Fwin/2; 0291 fhi=Pxx_mean_PeakFreq+Fwin/2; 0292 sel2=(flo<=f & f<=fhi); 0293 0294 % calculate area under spectrum around peak 0295 Pxx_mean_PeakArea = sum(tmpPxx(sel2))*(f(2)-f(1)); 0296 else 0297 Pxx_mean_PeakFreq=nan; 0298 Pxx_mean_PeakArea=nan; 0299 end 0300 Pxx_mean=tmpPxx; 0301 end 0302 0303 %% Add resulting power spectra to data structure 0304 % organization scheme: 0305 % data.VARIABLE_Power_SUA.(Pxx,PeakFreq,PeakArea,frequency) 0306 % data.VARIABLE_Power_MUA.(Pxx,PeakFreq,PeakArea,frequency) 0307 data.([var '_Power_SUA' options.output_suffix]).Pxx=Pxx; 0308 data.([var '_Power_SUA' options.output_suffix]).PeakFreq=PeakFreq; 0309 data.([var '_Power_SUA' options.output_suffix]).PeakArea=PeakArea; 0310 data.([var '_Power_SUA' options.output_suffix]).frequency=f; 0311 data.([var '_Power_MUA' options.output_suffix]).Pxx=Pxx_mean; 0312 data.([var '_Power_MUA' options.output_suffix]).PeakFreq=Pxx_mean_PeakFreq; 0313 data.([var '_Power_MUA' options.output_suffix]).PeakArea=Pxx_mean_PeakArea; 0314 data.([var '_Power_MUA' options.output_suffix]).frequency=f; 0315 0316 if ~ismember([var '_Power_SUA' options.output_suffix],data.results) 0317 data.results{end+1}=[var '_Power_SUA' options.output_suffix]; 0318 end 0319 0320 if ~ismember([var '_Power_MUA' options.output_suffix],data.results) 0321 data.results{end+1}=[var '_Power_MUA' options.output_suffix]; 0322 end 0323 0324 if options.exclude_data_flag 0325 for l=1:length(data.labels) 0326 data=rmfield(data,data.labels{l}); 0327 end 0328 end 0329 % % alternate organization scheme: 0330 % data.([var '_Pxx'])=Pxx; 0331 % data.([var '_Pxx_PeakFreq'])=PeakFreq; 0332 % data.([var '_Pxx_PeakArea'])=PeakArea; 0333 % data.([var '_Pxx_mean'])=Pxx_mean; 0334 % data.([var '_Pxx_mean_PeakFreq'])=Pxx_mean_PeakFreq; 0335 % data.([var '_Pxx_mean_PeakArea'])=Pxx_mean_PeakArea; 0336 % if ~ismember([var '_Pxx'],data.results) 0337 % data.results{end+1}=[var '_Pxx']; 0338 % data.results{end+1}=[var '_Pxx_PeakFreq']; 0339 % data.results{end+1}=[var '_Pxx_PeakArea']; 0340 % data.results{end+1}=[var '_Pxx_mean']; 0341 % data.results{end+1}=[var '_Pxx_mean_Pxx_PeakFreq']; 0342 % data.results{end+1}=[var '_Pxx_mean_Pxx_PeakArea']; 0343 % end 0344 0345 end 0346 0347 %% auto_gen_test_data_flag argout 0348 if options.auto_gen_test_data_flag 0349 argout = {data}; % specific to this function 0350 0351 dsUnitSaveAutoGenTestData(argin, argout); 0352 end