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