1 |
function [SM,EP]=dirspec(ID,SM,EP,varargin); |
---|
2 |
|
---|
3 |
%DIWASP V1.1 function with GD edit |
---|
4 |
%dirspec: main spectral estimation routine |
---|
5 |
% |
---|
6 |
%Note that there have been some chages made to this function since V1.1 |
---|
7 |
% |
---|
8 |
%-The outdated matlab function csd was replaced with the new version cpsd on |
---|
9 |
%line 90 |
---|
10 |
% |
---|
11 |
%-The ability to use along beam radial velocities has been added, and so the transfer |
---|
12 |
%parameters have been edited (line 109) so that the x,y position of the data input |
---|
13 |
%is also included |
---|
14 |
% |
---|
15 |
%-In addition, the first data input cell is radial velocity data, the initial 1-d |
---|
16 |
%frequency spectrum will be generated using the spectral average of all of |
---|
17 |
%the data input. Note that if the first data input cell is something other |
---|
18 |
%than 'radial' (pressure or sea surface height may be preferable), |
---|
19 |
%the initial 1-d frequency spectrum will be generated using ONLY that first |
---|
20 |
%data input. (see line 123) |
---|
21 |
% |
---|
22 |
% |
---|
23 |
%[SMout,EPout]=dirspec(ID,SM,EP,{options}) |
---|
24 |
% |
---|
25 |
%Outputs: |
---|
26 |
%SMout A spectral matrix structure containing the results |
---|
27 |
%Epout The estimation parameters structure with the values actually used for the computation including any default settings. |
---|
28 |
% |
---|
29 |
%Inputs: |
---|
30 |
%ID An instrument data structure containing the measured data |
---|
31 |
%SM A spectral matrix structure; data in field SM.S is ignored. |
---|
32 |
%EP The estimation parameters structure. For default values enter EP as [] |
---|
33 |
%[options] options entered as cell array with parameter/value pairs: e.g.{'MESSAGE',1,'PLOTTYPE',2}; |
---|
34 |
% Available options with default values: |
---|
35 |
% 'MESSAGE',1, Level of screen display: 0,1,2 (increasing output) |
---|
36 |
% 'PLOTTYPE',1, Plot type: 0 none, 1 3d surface, 2 polar type plot, 3 3d surface(compass angles), 4 polar plot(compass angles) |
---|
37 |
% 'FILEOUT','' Filename for output file: empty string means no file output |
---|
38 |
% |
---|
39 |
%Input structures ID and SM are required. Either [EP] or [options] can be included but must be in order if both are included. |
---|
40 |
%Type:%"help data_structures" for information on the DIWASP data structures |
---|
41 |
|
---|
42 |
%All of the implemented calculation algorithms are as described by: |
---|
43 |
%Hashimoto,N. 1997 "Analysis of the directional wave spectrum from field data" |
---|
44 |
%In: Advances in Coastal Engineering Vol.3. Ed:Liu,P.L-F. Pub:World Scientific,Singapore |
---|
45 |
% |
---|
46 |
%Copyright (C) 2002 Coastal Oceanography Group, CWR, UWA, Perth |
---|
47 |
|
---|
48 |
Options=struct('MESSAGE',1,'PLOTTYPE',1,'FILEOUT',''); |
---|
49 |
|
---|
50 |
if nargin<3 |
---|
51 |
error('All inputs other than OPTIONS required'); |
---|
52 |
elseif nargin>=4 |
---|
53 |
nopts=length(varargin{1}); |
---|
54 |
end |
---|
55 |
|
---|
56 |
ID=check_data(ID,1);if isempty(ID) return;end; |
---|
57 |
SM=check_data(SM,2);if isempty(SM) return;end; |
---|
58 |
EP=check_data(EP,3);if isempty(EP) return;end; |
---|
59 |
|
---|
60 |
if ~isempty(nopts) |
---|
61 |
if(rem(nopts,2)~=0) |
---|
62 |
warning('Options must be in Name/Value pairs - setting to defaults'); |
---|
63 |
else |
---|
64 |
for i=1:(nopts/2) |
---|
65 |
arg=varargin{1}{(2*i)}; |
---|
66 |
field=varargin{1}{(2*i-1)}; |
---|
67 |
Options=setfield(Options,field,arg);end;end |
---|
68 |
end |
---|
69 |
|
---|
70 |
ptype=Options.PLOTTYPE;displ=Options.MESSAGE; |
---|
71 |
|
---|
72 |
disp(' ');disp('calculating.....');disp(' ');disp('cross power spectra'); |
---|
73 |
|
---|
74 |
data=detrend(ID.data); |
---|
75 |
szd=size(ID.data,2); |
---|
76 |
dt=1/(ID.fs); |
---|
77 |
|
---|
78 |
%get resolution of FFT - if not specified, calculate a sensible value depending on sampling frequency |
---|
79 |
if isempty(EP.nfft) |
---|
80 |
nfft=2^(8+round(log2(ID.fs))); |
---|
81 |
EP.nfft=nfft; |
---|
82 |
else |
---|
83 |
nfft=EP.nfft; |
---|
84 |
end |
---|
85 |
|
---|
86 |
% NOTE: changed from csd to cpsd 6/15/07 |
---|
87 |
%calculate the cross-power spectra |
---|
88 |
for m=1:szd |
---|
89 |
for n=1:szd |
---|
90 |
[xpstmp,Ftmp]=cpsd(data(:,m),data(:,n),[],[],nfft,ID.fs); |
---|
91 |
xps(m,n,:)=xpstmp(2:(nfft/2)+1); |
---|
92 |
end |
---|
93 |
end |
---|
94 |
|
---|
95 |
F=Ftmp(2:(nfft/2)+1);nf=nfft/2; |
---|
96 |
|
---|
97 |
%calculate wavenumbers |
---|
98 |
disp('wavenumbers') |
---|
99 |
wns=wavenumber(2*pi*F,ID.depth*ones(size(F))); |
---|
100 |
pidirs=[-pi:(2*pi/EP.dres):pi-(2*pi/EP.dres)]; % this is dirs |
---|
101 |
|
---|
102 |
% NOTE: edited the transfer parameters to include x and y position |
---|
103 |
% x is (layout(1,m)),y is (layout(2,m)) (6/15/07) |
---|
104 |
|
---|
105 |
%calculate transfer parameters |
---|
106 |
disp('transfer parameters'); |
---|
107 |
disp(' '); |
---|
108 |
for m=1:szd |
---|
109 |
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 |
---|
110 |
for n=1:szd |
---|
111 |
kx(m,n,:,:)=wns*((ID.layout(1,n)-ID.layout(1,m))*cos(pidirs)+(ID.layout(2,n)-ID.layout(2,m))*sin(pidirs)); |
---|
112 |
end |
---|
113 |
end |
---|
114 |
|
---|
115 |
for m=1:szd |
---|
116 |
tfn(:,:)=trm(m,:,:); |
---|
117 |
Sxps(1:nf)=2*dt*xps(m,m,:); |
---|
118 |
Ss(m,:)=Sxps./(max(tfn').*conj(max(tfn'))); |
---|
119 |
end |
---|
120 |
|
---|
121 |
%have Ss = the spectral average if the radial velocities are used without a |
---|
122 |
%pressure or vertical range datatype first |
---|
123 |
compare=strcmp('radial',ID.datatypes(1)); |
---|
124 |
|
---|
125 |
if compare == 1 |
---|
126 |
Ss=mean(Ss); |
---|
127 |
end |
---|
128 |
|
---|
129 |
ffs=sum(F<=max(SM.freqs)); |
---|
130 |
|
---|
131 |
% call appropriate estimation function |
---|
132 |
disp(['directional spectra using' blanks(1) EP.method ' method']);disp(' ') |
---|
133 |
Specout=feval(EP.method,xps(:,:,1:ffs),trm(:,1:ffs,:),kx(:,:,1:ffs,:),Ss(:,1:ffs),pidirs,EP.iter,displ); |
---|
134 |
|
---|
135 |
Specout(find(isnan(Specout)))=0.0; |
---|
136 |
Hs=Hsig(Specout,F(2)-F(1),pidirs(2)-pidirs(1)); |
---|
137 |
|
---|
138 |
% map spectrum onto user defined spectrum matrix - need extra line of frequencies to avoid NaNs |
---|
139 |
[df,ddir]=meshgrid(SM.freqs,SM.dirs); |
---|
140 |
pidirs(EP.dres+1)=pi; |
---|
141 |
Specout=Specout'; |
---|
142 |
Specout(EP.dres+1,:)=Specout(1,:); |
---|
143 |
[Ff,Dd]=meshgrid(F(1:ffs),(180/pi)*pidirs); |
---|
144 |
S=interp2(Ff,Dd,Specout,df,ddir,'nearest'); |
---|
145 |
S=S*pi/180; |
---|
146 |
S(find(isnan(S)))=0.0; |
---|
147 |
|
---|
148 |
%check Hsig of mapped spectrum and check sufficiently close to original |
---|
149 |
Hs2=Hsig(S,SM.freqs(2)-SM.freqs(1),SM.dirs(2)-SM.dirs(1)); |
---|
150 |
if (Hs2-Hs)/Hs >0.01 |
---|
151 |
warning('User defined grid may be too coarse; try increasing resolution of ''SM.freqs'' or ''SM.dirs'''); |
---|
152 |
end |
---|
153 |
|
---|
154 |
%smooth spectrum |
---|
155 |
if(strcmp(EP.smooth,'ON')) |
---|
156 |
disp(' ');disp('smoothing spectrum...');disp(' '); |
---|
157 |
S=smoothspec(S,[1 0.5 0.25;1 0.5 0.25]); |
---|
158 |
end |
---|
159 |
SM.S=S'; |
---|
160 |
|
---|
161 |
infospec(SM); |
---|
162 |
|
---|
163 |
%write out spectrum matrix in DIWASP format |
---|
164 |
filename=Options.FILEOUT; |
---|
165 |
if(size(filename,2)>0) |
---|
166 |
disp('writing out spectrum matrix to file'); |
---|
167 |
writespec(SM,filename); |
---|
168 |
end |
---|
169 |
|
---|
170 |
%plot spectrum |
---|
171 |
if(ptype>0) |
---|
172 |
disp('finished...plotting spectrum'); |
---|
173 |
plotspec(SM,ptype); |
---|
174 |
T=['Directional spectrum estimated using ' blanks(1) EP.method ' method'];title(T); |
---|
175 |
end |
---|