Home > chronux_1_50 > 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: 
 Note units have to be consistent. See chronux.m for more information.
       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 form [NW K] e.g [3 5]) -- optional. If not 
                                                 specified, use [NW K]=[3 5]
            pad            (padding factor for the FFT) - optional (can take values -1,0,1,2...). 
                    -1 corresponds to no padding, 0 corresponds to padding
                    to the next highest power of 2 etc.
                       e.g. For N = 500, if PAD = -1, we do not pad; if PAD = 0, we pad the FFT
                       to 512 points, if pad=1, we pad to 1024 points 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: 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 form [NW K] e.g [3 5]) -- optional. If not
0027 %                                                 specified, use [NW K]=[3 5]
0028 %            pad            (padding factor for the FFT) - optional (can take values -1,0,1,2...).
0029 %                    -1 corresponds to no padding, 0 corresponds to padding
0030 %                    to the next highest power of 2 etc.
0031 %                       e.g. For N = 500, if PAD = -1, we do not pad; if PAD = 0, we pad the FFT
0032 %                       to 512 points, if pad=1, we pad to 1024 points etc.
0033 %                       Defaults to 0.
0034 %           Fs   (sampling frequency) - optional. Default 1.
0035 %           fpass    (frequency band to be used in the calculation in the form
0036 %                                   [fmin fmax])- optional.
0037 %                                   Default all frequencies between 0 and
0038 %                                   Fs/2
0039 % Output:
0040 %       SS       (second spectrum in form frequency x segments x trials x channels
0041 %                if segave=0; in the form frequency x trials x channels if segave=1)
0042 %       tau      (frequencies)
0043 if nargin < 3; error('Need data,segment duration and taper information'); end;
0044 if nargin < 4 ; params=[]; end;
0045 [tapers,pad,Fs,fpass,err,trialave,params]=getparams(params); 
0046 [N,Ntr,NC]=size(data);
0047 if Ntr==1; error('cannot compute second spectrum with just one trial'); end;
0048 dt=1/Fs; % sampling interval
0049 T=N*dt; % length of data in seconds
0050 E=0:win:T-win; % fictitious event triggers
0051 datatmp=createdatamatc(data(:,1,1),E,Fs,[0 win]); % segmented data
0052 Ninseg=size(datatmp,1); % number of samples in segments
0053 nfft=max(2^(nextpow2(Ninseg)+pad),Ninseg);
0054 [f,findx]=getfgrid(Fs,nfft,fpass); 
0055 NF=length(findx);
0056 S=zeros(NF,Ntr,NC);
0057 for nc=1:NC;
0058     for ntr=1:Ntr;
0059         datatmp=change_row_to_column(data(:,ntr,nc));
0060         s=mtspectrumsegc(datatmp,win,params,1);
0061         S(:,ntr,nc)=s;
0062     end
0063 end;
0064 Sm=mean(S,2);
0065 if use==1;
0066    params.tapers=tapers2spec;
0067    params.Fs=1/(f(2)-f(1));
0068    params.fpass=[0 params.Fs/2];
0069 else;
0070    tau=[0:NF-1]/max(f);
0071    cosinefunc=cos(2*pi*f'*tau);
0072 end;
0073 
0074 for nc=1:NC;
0075     for ntr=1:Ntr;
0076         s=S(:,ntr,nc)./Sm(:,nc);
0077         s=log(s);
0078         if use==1;
0079             sflip=flipdim(s,1);
0080             s=[sflip(1:NF-1);s];
0081             [ss,tau]=mtspectrumc(s,params);
0082             SS(:,ntr,nc)=ss;
0083         else;
0084             s=repmat(s,[1 NF]).*cosinefunc;
0085     %         subplot(221); plot(s(:,1));
0086     %         subplot(222); plot(s(:,10));
0087     %         subplot(223); plot(s(:,100));
0088     %         subplot(224); plot(s(:,120));
0089     %         pause
0090             s=trapz(f,s,1)';
0091             ss=s.*conj(s);
0092 %         plot(tau,s)
0093 %         pause
0094         end
0095         SS(:,ntr,nc)=ss;
0096     end
0097 end;
0098 SS=mean(SS,2);

Generated on Mon 09-Oct-2006 00:54:52 by m2html © 2003