MTSpectrum - Compute LFP spectrum by multi-taper estimation.
The spectrum is computed as the average spectrogram, using the <a href="http://www.chronux.org">chronux</a> toolbox.
USAGE
[spectrum,f,s] = MTSpectrum(lfp,<options>)
lfp unfiltered LFP <a href="matlab:help samples">samples</a> (one channel).
<options> optional list of property-value pairs (see table below)
=========================================================================
Properties Values
-------------------------------------------------------------------------
'frequency' sampling rate (in Hz) (default = from timestamps if
available, otherwise 1250Hz)
'range' frequency range (in Hz) (default = all)
'window' duration (in s) of the time window (default = 5)
'overlap' overlap between successive windows (default = window/2)
'step' step between successive windows (default = window/2)
'tapers' relative resolution and order of the tapers [NW K]
(default = [3 5])
'pad' FFT padding (see help for <a href="matlab:help mtspecgramc">mtspecgramc</a>) (default = 0)
'show' plot log spectrum (default = 'off')
=========================================================================
NOTES
The LFP can be provided either as a time stamped matrix (list of time-voltage
pairs), or as a voltage vector - in which case the frequency must be specified.
The time displacement between successive short time spectra can be supplied
either as a 'step' (explicit time difference) or as an 'overlap' (between
successive time windows).
OUTPUT
spectrum spectrum
f frequency bins
s error
DEPENDENCIES
This function requires the <a href="http://www.chronux.org">chronux</a> toolbox.
SEE
See also MTSpectrogram, SpectrogramBands, MTCoherence, MTCoherogram, PlotMean.0001 function [spectrum,f,s] = MTSpectrum(lfp,varargin) 0002 0003 %MTSpectrum - Compute LFP spectrum by multi-taper estimation. 0004 % 0005 % The spectrum is computed as the average spectrogram, using the <a href="http://www.chronux.org">chronux</a> toolbox. 0006 % 0007 % USAGE 0008 % 0009 % [spectrum,f,s] = MTSpectrum(lfp,<options>) 0010 % 0011 % lfp unfiltered LFP <a href="matlab:help samples">samples</a> (one channel). 0012 % <options> optional list of property-value pairs (see table below) 0013 % 0014 % ========================================================================= 0015 % Properties Values 0016 % ------------------------------------------------------------------------- 0017 % 'frequency' sampling rate (in Hz) (default = from timestamps if 0018 % available, otherwise 1250Hz) 0019 % 'range' frequency range (in Hz) (default = all) 0020 % 'window' duration (in s) of the time window (default = 5) 0021 % 'overlap' overlap between successive windows (default = window/2) 0022 % 'step' step between successive windows (default = window/2) 0023 % 'tapers' relative resolution and order of the tapers [NW K] 0024 % (default = [3 5]) 0025 % 'pad' FFT padding (see help for <a href="matlab:help mtspecgramc">mtspecgramc</a>) (default = 0) 0026 % 'show' plot log spectrum (default = 'off') 0027 % ========================================================================= 0028 % 0029 % NOTES 0030 % 0031 % The LFP can be provided either as a time stamped matrix (list of time-voltage 0032 % pairs), or as a voltage vector - in which case the frequency must be specified. 0033 % 0034 % The time displacement between successive short time spectra can be supplied 0035 % either as a 'step' (explicit time difference) or as an 'overlap' (between 0036 % successive time windows). 0037 % 0038 % OUTPUT 0039 % 0040 % spectrum spectrum 0041 % f frequency bins 0042 % s error 0043 % 0044 % DEPENDENCIES 0045 % 0046 % This function requires the <a href="http://www.chronux.org">chronux</a> toolbox. 0047 % 0048 % SEE 0049 % 0050 % See also MTSpectrogram, SpectrogramBands, MTCoherence, MTCoherogram, PlotMean. 0051 0052 % Copyright (C) 2010-2014 by Michaël Zugaro 0053 % 0054 % This program is free software; you can redistribute it and/or modify 0055 % it under the terms of the GNU General Public License as published by 0056 % the Free Software Foundation; either version 3 of the License, or 0057 % (at your option) any later version. 0058 0059 % Make sure chronux is installed and functional 0060 CheckChronux('mtspecgramc'); 0061 0062 % Defaults 0063 f = 1250; 0064 frequency = []; 0065 window = 5; 0066 range = []; 0067 overlap = []; 0068 step = []; 0069 show = 'off'; 0070 tapers = [3 5]; 0071 pad = 0; 0072 0073 % Check number of parameters 0074 if nargin < 1 | mod(length(varargin),2) ~= 0, 0075 error('Incorrect number of parameters (type ''help <a href="matlab:help MTSpectrum">MTSpectrum</a>'' for details).'); 0076 end 0077 0078 % Check parameter sizes 0079 if size(lfp,2) ~= 1 && size(lfp,2) ~= 2, 0080 error('Parameter ''lfp'' is not a vector or a Nx2 matrix (type ''help <a href="matlab:help MTSpectrum">MTSpectrum</a>'' for details).'); 0081 end 0082 0083 % Parse parameter list 0084 v = {}; 0085 for i = 1:2:length(varargin), 0086 if ~ischar(varargin{i}), 0087 error(['Parameter ' num2str(i+2) ' is not a property (type ''help <a href="matlab:help MTSpectrum">MTSpectrum</a>'' for details).']); 0088 end 0089 switch(lower(varargin{i})), 0090 case 'frequency', 0091 frequency = varargin{i+1}; 0092 if ~isdscalar(frequency,'>0'), 0093 error('Incorrect value for property ''frequency'' (type ''help <a href="matlab:help MTSpectrum">MTSpectrum</a>'' for details).'); 0094 end 0095 case 'range', 0096 range = varargin{i+1}; 0097 if ~isdvector(range,'#2','<','>=0'), 0098 error('Incorrect value for property ''range'' (type ''help <a href="matlab:help MTSpectrum">MTSpectrum</a>'' for details).'); 0099 end 0100 case 'window', 0101 window = varargin{i+1}; 0102 if ~isdscalar(window,'>0'), 0103 error('Incorrect value for property ''window'' (type ''help <a href="matlab:help MTSpectrum">MTSpectrum</a>'' for details).'); 0104 end 0105 case 'overlap', 0106 overlap = varargin{i+1}; 0107 if ~isdscalar(overlap,'>0'), 0108 error('Incorrect value for property ''overlap'' (type ''help <a href="matlab:help MTSpectrum">MTSpectrum</a>'' for details).'); 0109 end 0110 case 'step', 0111 step = varargin{i+1}; 0112 if ~isdscalar(step,'>0'), 0113 error('Incorrect value for property ''step'' (type ''help <a href="matlab:help MTSpectrum">MTSpectrum</a>'' for details).'); 0114 end 0115 case 'tapers', 0116 tapers = varargin{i+1}; 0117 if ~isivector(tapers,'#2','>0'), 0118 error('Incorrect value for property ''tapers'' (type ''help <a href="matlab:help MTSpectrum">MTSpectrum</a>'' for details).'); 0119 end 0120 case 'pad', 0121 pad = varargin{i+1}; 0122 if ~isiscalar(pad,'>-1'), 0123 error('Incorrect value for property ''pad'' (type ''help <a href="matlab:help MTSpectrum">MTSpectrum</a>'' for details).'); 0124 end 0125 case 'show', 0126 show = varargin{i+1}; 0127 if ~isastring(show,'on','off'), 0128 error('Incorrect value for property ''show'' (type ''help <a href="matlab:help MTSpectrum">MTSpectrum</a>'' for details).'); 0129 end 0130 otherwise, 0131 error(['Unknown property ''' num2str(varargin{i}) ''' (type ''help <a href="matlab:help MTSpectrum">MTSpectrum</a>'' for details).']); 0132 end 0133 if ~strcmp(varargin{i},'show'), v = {v{:},varargin{i:i+1}}; end 0134 end 0135 0136 % Compute spectrogram and moments 0137 [spectrogram,~,f] = MTSpectrogram(lfp,v{:}); 0138 spectrogram = spectrogram'; 0139 mu = mean(spectrogram); 0140 v = var(spectrogram); 0141 0142 % Compute spectrum 0143 spectrum = mu; 0144 s = sqrt(v); 0145 0146 % Plot log, i.e. mean E[log(spectrogram)] and stdev sqrt(Var[log(spectrogram)]) 0147 % (see http://en.wikipedia.org/wiki/Taylor_expansions_for_the_moments_of_functions_of_random_variables) 0148 if strcmp(lower(show),'on'), 0149 figure; 0150 %logSpectrum = log(spectrum); % classic value 0151 logSpectrum = log(mu)-v./(2*mu.*mu); % corrected value, valid for mu/s > 1.5 0152 logS = sqrt(v./mu.^2); % valid for mu/s > 2.5 0153 PlotMean(f,logSpectrum,logSpectrum-logS,logSpectrum+logS,':'); 0154 xlabel('Frequency (Hz)'); 0155 ylabel('Power'); 0156 title('Power Spectrum'); 0157 end