# 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 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...).
                        -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:

- 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
- change_row_to_column Helper routine to transform 1d arrays into column vectors that are needed
- getfgrid Helper function that gets the frequency grid associated with a given fft based computation
- getparams Helper function to convert structure params to variables used by the

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...).
    0045 %                    -1 corresponds to no padding, 0 corresponds to padding
    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;
    0061 [tapers,pad,Fs,fpass,err,trialave,params]=getparams(params);
    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
    0069 nfft=max(2^(nextpow2(Ninseg)+pad),Ninseg);
    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);