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