Home > chronux > spectral_analysis > continuous > mtspectrum_of_spectrumc.m

# mtspectrum_of_spectrumc

## PURPOSE Multi-taper segmented, second spectrum (spectrum of the log spectrum) for a continuous process

## SYNOPSIS function [SS,tau]=mtspectrum_of_spectrumc(data,win,tapers2spec,params)

## DESCRIPTION ``` Multi-taper segmented, second spectrum (spectrum of the log spectrum) for a continuous process
This routine computes the second spectrum by explicitly evaluating the
Fourier transform (since the spectrum is symmetric in frequency, it uses
a cosine transform)

Usage:

[SS,tau]=mtspectrum_of_spectrumc(data,win,tapers2spec,params)
Input:
data (single channel) -- required
win  (duration of the segments) - required.
tapers2spec (tapers used for the spectrum of spectrum computation) -
required in the form [use TW K] - Note that spectrum of the
spectrum involves computing two Fourier transforms. While the first
transform (of the original data) is always computed using the
multi-taper method, the current routine allows the user to specify
whether or not to use this method for the second transform. use=1
means use tapers, use=anything other than 1 means do not use the
multitaper method. If use=1, then tapers2spec controls the
smoothing for the second Fourier transform. Otherwise, a direct
Fourier transform is computed.
params: structure with fields tapers, pad, Fs, fpass, err, trialave
- optional
tapers : precalculated tapers from dpss or in the one of the following
forms:
(1) A numeric vector [TW K] where TW is the
time-bandwidth product and K is the number of
tapers to be used (less than or equal to
2TW-1).
(2) A numeric vector [W T p] where W is the
bandwidth, T is the duration of the data and p
is an integer such that 2TW-p tapers are used. In
this form there is no default i.e. to specify
the bandwidth, you have to specify T and p as
well. Note that the units of W and T have to be
consistent: if W is in Hz, T must be in seconds
and vice versa. Note that these units must also
be consistent with the units of params.Fs: W can
be in Hz if and only if params.Fs is in Hz.
The default is to use form 1 with TW=3 and K=5

pad            (padding factor for the FFT) - optional (can take values -1,0,1,2...).
to the next highest power of 2 etc.
Defaults to 0.
Fs   (sampling frequency) - optional. Default 1.
fpass    (frequency band to be used in the calculation in the form
[fmin fmax])- optional.
Default all frequencies between 0 and
Fs/2
Output:
SS       (second spectrum in form frequency x segments x trials x channels
if segave=0; in the form frequency x trials x channels if segave=1)
tau      (frequencies)```

## CROSS-REFERENCE INFORMATION This function calls:
• createdatamatc Helper function to create an event triggered matrix from univariate
• mtspectrumc Multi-taper spectrum - continuous process
• mtspectrumsegc Multi-taper segmented spectrum for a univariate continuous process
This function is called by:

## SOURCE CODE ```0001 function [SS,tau]=mtspectrum_of_spectrumc(data,win,tapers2spec,params)
0002 % Multi-taper segmented, second spectrum (spectrum of the log spectrum) for a continuous process
0003 % This routine computes the second spectrum by explicitly evaluating the
0004 % Fourier transform (since the spectrum is symmetric in frequency, it uses
0005 % a cosine transform)
0006 %
0007 % Usage:
0008 %
0009 % [SS,tau]=mtspectrum_of_spectrumc(data,win,tapers2spec,params)
0010 % Input:
0011 % Note units have to be consistent. See chronux.m for more information.
0012 %       data (single channel) -- required
0013 %       win  (duration of the segments) - required.
0014 %       tapers2spec (tapers used for the spectrum of spectrum computation) -
0015 %       required in the form [use TW K] - Note that spectrum of the
0016 %       spectrum involves computing two Fourier transforms. While the first
0017 %       transform (of the original data) is always computed using the
0018 %       multi-taper method, the current routine allows the user to specify
0019 %       whether or not to use this method for the second transform. use=1
0020 %       means use tapers, use=anything other than 1 means do not use the
0021 %       multitaper method. If use=1, then tapers2spec controls the
0022 %       smoothing for the second Fourier transform. Otherwise, a direct
0023 %       Fourier transform is computed.
0024 %       params: structure with fields tapers, pad, Fs, fpass, err, trialave
0025 %       - optional
0026 %           tapers : precalculated tapers from dpss or in the one of the following
0027 %                    forms:
0028 %                    (1) A numeric vector [TW K] where TW is the
0029 %                        time-bandwidth product and K is the number of
0030 %                        tapers to be used (less than or equal to
0031 %                        2TW-1).
0032 %                    (2) A numeric vector [W T p] where W is the
0033 %                        bandwidth, T is the duration of the data and p
0034 %                        is an integer such that 2TW-p tapers are used. In
0035 %                        this form there is no default i.e. to specify
0036 %                        the bandwidth, you have to specify T and p as
0037 %                        well. Note that the units of W and T have to be
0038 %                        consistent: if W is in Hz, T must be in seconds
0039 %                        and vice versa. Note that these units must also
0040 %                        be consistent with the units of params.Fs: W can
0041 %                        be in Hz if and only if params.Fs is in Hz.
0042 %                        The default is to use form 1 with TW=3 and K=5
0043 %
0044 %            pad            (padding factor for the FFT) - optional (can take values -1,0,1,2...).
0046 %                    to the next highest power of 2 etc.
0047 %                       e.g. For N = 500, if PAD = -1, we do not pad; if PAD = 0, we pad the FFT
0048 %                       to 512 points, if pad=1, we pad to 1024 points etc.
0049 %                       Defaults to 0.
0050 %           Fs   (sampling frequency) - optional. Default 1.
0051 %           fpass    (frequency band to be used in the calculation in the form
0052 %                                   [fmin fmax])- optional.
0053 %                                   Default all frequencies between 0 and
0054 %                                   Fs/2
0055 % Output:
0056 %       SS       (second spectrum in form frequency x segments x trials x channels
0057 %                if segave=0; in the form frequency x trials x channels if segave=1)
0058 %       tau      (frequencies)
0059 if nargin < 3; error('Need data,segment duration and taper information'); end;
0060 if nargin < 4 ; params=[]; end;
0062 [N,Ntr,NC]=size(data);
0063 if Ntr==1; error('cannot compute second spectrum with just one trial'); end;
0064 dt=1/Fs; % sampling interval
0065 T=N*dt; % length of data in seconds
0066 E=0:win:T-win; % fictitious event triggers
0067 datatmp=createdatamatc(data(:,1,1),E,Fs,[0 win]); % segmented data
0068 Ninseg=size(datatmp,1); % number of samples in segments
0070 [f,findx]=getfgrid(Fs,nfft,fpass);
0071 NF=length(findx);
0072 S=zeros(NF,Ntr,NC);
0073 for nc=1:NC;
0074     for ntr=1:Ntr;
0075         datatmp=change_row_to_column(data(:,ntr,nc));
0076         s=mtspectrumsegc(datatmp,win,params,1);
0077         S(:,ntr,nc)=s;
0078     end
0079 end;
0080 Sm=mean(S,2);
0081 if use==1;
0082    params.tapers=tapers2spec;
0083    params.Fs=1/(f(2)-f(1));
0084    params.fpass=[0 params.Fs/2];
0085 else;
0086    tau=[0:NF-1]/max(f);
0087    cosinefunc=cos(2*pi*f'*tau);
0088 end;
0089
0090 for nc=1:NC;
0091     for ntr=1:Ntr;
0092         s=S(:,ntr,nc)./Sm(:,nc);
0093         s=log(s);
0094         if use==1;
0095             sflip=flipdim(s,1);
0096             s=[sflip(1:NF-1);s];
0097             [ss,tau]=mtspectrumc(s,params);
0098             SS(:,ntr,nc)=ss;
0099         else;
0100             s=repmat(s,[1 NF]).*cosinefunc;
0101     %         subplot(221); plot(s(:,1));
0102     %         subplot(222); plot(s(:,10));
0103     %         subplot(223); plot(s(:,100));
0104     %         subplot(224); plot(s(:,120));
0105     %         pause
0106             s=trapz(f,s,1)';
0107             ss=s.*conj(s);
0108 %         plot(tau,s)
0109 %         pause
0110         end
0111         SS(:,ntr,nc)=ss;
0112     end
0113 end;
0114 SS=mean(SS,2);```

Generated on Fri 28-Sep-2012 12:34:30 by m2html © 2005