


Space frequency SVD of input data - continuous processes
Usage: [sv,sp,fm] = spsvd(data,tapers,Fs,fpass,mdkp)
Inputs:
data (data matrix in timexchannels form)-required
tapers (tapers parameters or precalculated tapers)-required
Fs (sampling frequency)-optional. Default 1
fpass (band of frequencies to be kept)-optional. Default [0 Fs/2]
pad (padding factor for the FFT) - optional. Defaults to 0.
e.g. For N = 500, if PAD = 0, we pad the FFT
to 512 points; if PAD = 2, we pad the FFT
to 2048 points, etc.
mdkp (number of dimensions to be kept)-optional. Default is the
maximum possible modes determined by taper parameters
Outputs:
sv sp fm : singular values, space modes, frequency modes

0001 function [sv,sp,fm] = spsvd(data,tapers,Fs,fpass,pad,mdkp) 0002 % Space frequency SVD of input data - continuous processes 0003 % Usage: [sv,sp,fm] = spsvd(data,tapers,Fs,fpass,mdkp) 0004 % Inputs: 0005 % data (data matrix in timexchannels form)-required 0006 % tapers (tapers parameters or precalculated tapers)-required 0007 % Fs (sampling frequency)-optional. Default 1 0008 % fpass (band of frequencies to be kept)-optional. Default [0 Fs/2] 0009 % pad (padding factor for the FFT) - optional. Defaults to 0. 0010 % e.g. For N = 500, if PAD = 0, we pad the FFT 0011 % to 512 points; if PAD = 2, we pad the FFT 0012 % to 2048 points, etc. 0013 % mdkp (number of dimensions to be kept)-optional. Default is the 0014 % maximum possible modes determined by taper parameters 0015 % 0016 % Outputs: 0017 % sv sp fm : singular values, space modes, frequency modes 0018 0019 0020 if nargin < 2; error('Need data and tapers'); end; 0021 if nargin < 3; Fs=1; end; 0022 if nargin < 4; fpass=[0 Fs/2]; end; 0023 0024 N=size(data,1); 0025 tapers=dpsschk(tapers,N); 0026 nfft=2^(nextpow2(N)+pad);% number of points in fft 0027 [N,K]=size(tapers); 0028 if isempty(mdkp) | nargin<5; mdkp=min(K,NCHAN); 0029 elseif mdkp > min(K,NCHAN); error('mdkp has to be less than both K and NCHAN');end; 0030 0031 tvec=1:N; 0032 tvec=tvec*2*pi*i; 0033 [f,findx]=getfgrid(Fs,nfft,fpass); 0034 nf=length(f); 0035 0036 sp=zeros(nwin,NCHAN,nf,mdkp); 0037 sp=sp+i*sp; 0038 fm=zeros(nwin,K,nf,mdkp); 0039 fm=fm+i*fm; 0040 sv=zeros(nwin,nf,min([K,NCHAN])); 0041 proj=zeros(N,K); 0042 for j=1:nf 0043 f0=f(j)/Fs; 0044 for k=1:K 0045 proj(:,k)=tapers(:,k).*exp(-f0*tvec'); 0046 end 0047 tmp=data'*proj; % projected data 0048 [u,s,v]= svd(tmp,0); % svd 0049 for mk=1:mdkp, 0050 sp(n,:,j,mk)=u(:,mk)'; 0051 fm(n,:,j,mk)=v(:,mk)'; 0052 end 0053 sv(n,j,:)=diag(s); 0054 end;