root/DPWP/trunk/DPWP/diwasp_1_1GD/dirspec.m

function [SM,EP]=dirspec(ID,SM,EP,varargin)

%DIWASP V1.1 function with GD edit
%dirspec: main spectral estimation routine
%
%IMPT!:  See code for recent changes as of 2/18/09
%
%Note that there have been some chages made to this function since V1.1
%
%-The outdated matlab function csd was replaced with the new version cpsd 
%
%
%-The ability to use along beam radial velocities has been added, and so the transfer 
%parameters have been edited so that the x,y position of the data input
%is also included
%
%-In addition, the first data input cell is radial velocity data, the initial 1-d
%frequency spectrum will be generated using the spectral average of all of
%the data input.  Note that if the first data input cell is something other
%than 'radial' (pressure or sea surface height may be preferable), 
%the initial 1-d frequency spectrum will be generated using ONLY that first
%data input. 
%
%-the way the cpsd function created the window length was changed to pick
%an option that better suits the data length used
%
%-The formulation of pidirs was changed from hardwired at -pi to pi, to bounding 
%by the user selected SM.dirs data field
%
%[SMout,EPout]=dirspec(ID,SM,EP,{options})
%
%Outputs:
%SMout      A spectral matrix structure containing the results
%Epout          The estimation parameters structure with the values actually used for the computation including any default settings.
%
%Inputs:
%ID                     An instrument data structure containing the measured data
%SM             A spectral matrix structure; data in field SM.S is ignored.
%EP                 The estimation parameters structure. For default values enter EP as []
%[options]  options entered as cell array with parameter/value pairs: e.g.{'MESSAGE',1,'PLOTTYPE',2};
%                Available options with default values:
%                    'MESSAGE',1,    Level of screen display: 0,1,2 (increasing output)
%                    'PLOTTYPE',1,   Plot type: 0 none, 1 3d surface, 2 polar type plot, 3 3d surface(compass angles), 4 polar plot(compass angles)
%                    'FILEOUT',''        Filename for output file: empty string means no file output
%                
%Input structures ID and SM are required. Either [EP] or [options] can be included but must be in order if both are included.
%Type:%"help data_structures" for information on the DIWASP data structures

%All of the implemented calculation algorithms are as described by:
%Hashimoto,N. 1997 "Analysis of the directional wave spectrum from field data" 
%In: Advances in Coastal Engineering Vol.3. Ed:Liu,P.L-F. Pub:World Scientific,Singapore
%
%Copyright (C) 2002 Coastal Oceanography Group, CWR, UWA, Perth

Options=struct('MESSAGE',1,'PLOTTYPE',1,'FILEOUT','');

if nargin<3
   error('All inputs other than OPTIONS required');
elseif nargin>=4
   nopts=length(varargin{1});
end

ID=check_data(ID,1);if isempty(ID) return;end;
SM=check_data(SM,2);if isempty(SM) return;end;
EP=check_data(EP,3);if isempty(EP) return;end;

if ~isempty(nopts)
if(rem(nopts,2)~=0)
        warning('Options must be in Name/Value pairs - setting to defaults');
        else
        for i=1:(nopts/2)
        arg=varargin{1}{(2*i)};
        field=varargin{1}{(2*i-1)};
        Options=setfield(Options,field,arg);
    end
end   
end

ptype=Options.PLOTTYPE;displ=Options.MESSAGE;

disp(' ');disp('calculating.....');disp(' ');disp('cross power spectra');

data=detrend(ID.data);
szd=size(ID.data,2);
dt=1/(ID.fs);

%get resolution of FFT - if not specified, calculate a sensible value depending on sampling frequency
if isempty(EP.nfft)
    nfft=2^(8+round(log2(ID.fs)));
    EP.nfft=nfft;
else
    nfft=EP.nfft;
end

%NOTE: changed on 9/5/08 window length for cpsd set to be a window equal to
%or slightly larger than nfft so padding is limited
dataLength=size(ID.data,1);
numWindows=ceil(dataLength/nfft);
window=ceil(dataLength/numWindows);

%Do a version check for the upcoming cpsd function. If the version is 2006b
%or later (date Aug,3,2006 or datenum 732892 cpsd will be handled
%differently.  added on 9/17/08
[v d]=version;
ver_date=datenum(d);


% NOTE: changed cross power spectral function from csd to cpsd 6/15/07
for m=1:szd
for n=1:szd
    % Need to have 2 different cpsd functions depending on Matlab version.
    % added of 9/17/08
    if ver_date < 732892
        [xpstmp,Ftmp]=cpsd(data(:,m),data(:,n),window,[],nfft,ID.fs);
    else
        [xpstmp,Ftmp]=cpsd(data(:,n),data(:,m),window,[],nfft,ID.fs);
    end;
   xps(m,n,:)=xpstmp(2:(nfft/2)+1);
   xpsout=xps;
end
end

F=Ftmp(2:(nfft/2)+1);nf=nfft/2;

%calculate wavenumbers
disp('wavenumbers')
wns=wavenumber(2*pi*F,ID.depth*ones(size(F)));
%Edited on 2/18/09 to adjust for whatever the user selects for the
%directional bounds (whatever SM.dirs is)
pidirs=[SM.dirs(1)*(pi/180):(2*pi/EP.dres):SM.dirs(end)*(pi/180)-(2*pi/EP.dres)]; 


%  NOTE:  edited the transfer parameters to include x and y position
%  x is (layout(1,m)),y is (layout(2,m))   (6/15/07)
%calculate transfer parameters
disp('transfer parameters');
disp(' ');
for m=1:szd
        [trm(m,:,:)]=feval(ID.datatypes{m},2*pi*F,pidirs,wns,ID.layout(3,m),ID.depth,ID.layout(1,m),ID.layout(2,m)); % 2piF is ffreqs
 for n=1:szd
    kx(m,n,:,:)=wns*((ID.layout(1,n)-ID.layout(1,m))*cos(pidirs)+(ID.layout(2,n)-ID.layout(2,m))*sin(pidirs));
end
end

for m=1:szd
   tfn(:,:)=trm(m,:,:);
   Sxps(1:nf)=2*dt*xps(m,m,:);
   Ss(m,:)=Sxps./(max(tfn').*conj(max(tfn')));
end
SsOut=Ss;

%What are the number of sensors (data series) used to compute the power
%spectrum?  If we aren't using radial velocities it is just 1.

dataUsed=1;


%NOTE: edited on 8/07
%have Ss = the spectral average if the radial velocities are used without a
%pressure or vertical range datatype first
lags=50;
compare=strcmp('radial',ID.datatypes(1));
if compare == 1
    Ss=mean(Ss);
    %since this is using radial velocities the number of data series used
    %the total number of sensors or szd
    dataUsed=szd;
    
    %also need to now calculate the effective number of windows to
    %accurately calculate the degrees of freedom below.
    Tau=ones(dataUsed); %Tau Matrix
    Nstar=Tau;  %Nstar Matrix
    for i=1:dataUsed
        for k=1:dataUsed
            covi=xcov(data(:,i),lags,'coeff');
            covk=xcov(data(:,k),lags,'coeff');
            covik=xcov(data(:,i),data(:,k),lags,'coeff');
            covki=xcov(data(:,k),data(:,i),lags,'coeff');
            covicovk=covi.*covk;
            covikcovki=covik.*covki;
            %sign1=sign(covicovk);
            %sign2=sign(covikcovki);
            %part1=(abs(covicovk)).^.5;
            %part2=(abs(covikcovki)).^.5;
            Tau(i,k)=sum(covicovk+covikcovki);
            Nstar=dataLength./Tau;
        end
    end
    %calculate the average value of Nstar for the off-diagnal values
    totalN=sum(sum(Nstar));
    for i=1:dataUsed
        out(i)=Nstar(i,i);
    end
    autoOut=sum(out); %don't count the auto values
    totalN=totalN-autoOut;
    avgN=totalN/((dataUsed^2)-dataUsed);%find the average N value for all co values
    Nstar=avgN;    
else
    %if there is only ONE time series used for the frequency spec
    covi=xcov(data(:,1),lags,'coeff');
    Tau=sum(covi.*covi);
    Nstar=dataLength./Tau;
end

%now find the effective number windows for each time series given Nstar
numWindows=Nstar/nfft;


%Compute the degrees of freedom depending on cpsd options and the data used
%See Emery and Thomson, 2004, added on 9/17/08
windowConstant=2.5164; %for Hamming Window (Matlab default for cpsd)
degF=2*numWindows*dataUsed*windowConstant;
SM.degF=degF;


ffs=sum(F<=max(SM.freqs));

    % call appropriate estimation function
disp(['directional spectra using' blanks(1) EP.method ' method']);disp(' ')

[Specout]=feval(EP.method,xps(:,:,1:ffs),trm(:,1:ffs,:),kx(:,:,1:ffs,:),Ss(:,1:ffs),pidirs,EP.iter,displ);

Specout(find(isnan(Specout)))=0.0;
Hs=Hsig(Specout,F(2)-F(1),pidirs(2)-pidirs(1));

% map spectrum onto user defined spectrum matrix - need extra line of frequencies to avoid NaNs
[df,ddir]=meshgrid(SM.freqs,SM.dirs);
pidirs(EP.dres+1)=SM.dirs(end); %edited on 2/18/09 to adjust for changes in pidirs (see above)
Specout=Specout';
Specout(EP.dres+1,:)=Specout(1,:);
[Ff,Dd]=meshgrid(F(1:ffs),(180/pi)*pidirs);
S=interp2(Ff,Dd,Specout,df,ddir,'nearest');
S=S*pi/180;
S(find(isnan(S)))=0.0;

%check Hsig of mapped spectrum and check sufficiently close to original
%*****As far as I can tell this warning is useless**** Edited out on 3/2/10

%Hs2=Hsig(S,SM.freqs(2)-SM.freqs(1),SM.dirs(2)-SM.dirs(1));
%if (Hs2-Hs)/Hs >0.01
%    warning('User defined grid may be too coarse; try increasing resolution of ''SM.freqs'' or ''SM.dirs''');
%end

%smooth spectrum
if(strcmp(EP.smooth,'ON'))
    disp(' ');disp('smoothing spectrum...');disp(' ');
    S=smoothspec(S,[1 0.5 0.25;1 0.5 0.25]);
end
SM.S=S';

infospec(SM);

%write out spectrum matrix in DIWASP format
filename=Options.FILEOUT;
if(size(filename,2)>0)
   disp('writing out spectrum matrix to file');
   writespec(SM,filename);
end

%plot spectrum
if(ptype>0)
    disp('finished...plotting spectrum');
    plotspec(SM,ptype);
    T=['Directional spectrum estimated using ' blanks(1) EP.method ' method'];title(T);
end