0001 function wspecplotwspecplot(S,varargin)
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023
0024
0025
0026
0027
0028
0029
0030
0031
0032
0033
0034
0035
0036
0037
0038
0039
0040
0041
0042
0043
0044
0045
0046
0047
0048
0049
0050
0051
0052
0053
0054
0055
0056
0057
0058
0059
0060
0061
0062
0063
0064
0065
0066
0067
0068
0069
0070
0071
0072
0073
0074
0075
0076
0077
0078
0079
0080
0081
0082
0083
0084
0085
0086
0087
0088
0089
0090
0091
0092
0093
0094
0095
0096
0097
0098 LegendOn = 1;
0099
0100 P = varargin;
0101
0102 ih=ishold;
0103 Ns=length(S);
0104 if Ns>1
0105 cfig=gcf;
0106 for ix=1:Ns,
0107 if ih
0108 newplot
0109 else
0110 figure(cfig-1+ix)
0111 end
0112 wspecplot(S(ix),P{:})
0113 end
0114 return
0115 end
0116
0117
0118
0119 plotflag = 1;
0120 lintype = 'b-';
0121 Nlevel = 7;
0122
0123 Np = length(P);
0124 if Np>0
0125 strix = zeros(1,Np);
0126 for ix=1:Np,
0127 strix(ix)=ischar(P{ix});
0128 end
0129 k = find(strix);
0130 if any(k)
0131 Nk = length(k);
0132 if Nk>1
0133 warning('More than 1 strings are not allowed')
0134 end
0135 lintype = P{k(1)};
0136 Np = Np-Nk;
0137 P = {P{find(~strix)}};
0138 end
0139 if Np>0 & ~isempty(P{1}), plotflag=P{1};end
0140 end
0141
0142
0143 ftype = freqtype(S);
0144 switch ftype
0145 case {'w'},
0146 freq = S.w;
0147 xlbl_txt = 'Frequency [rad/s]';
0148 ylbl1_txt = 'S(w) [m^2 s / rad]';
0149 ylbl3_txt = 'Directional Spectrum';
0150 zlbl_txt = 'S(w,\theta) [m^2 s / rad^2]';
0151 funit = ' [rad/s]';
0152 Sunit = ' [m^2 s / rad]';
0153 case {'f'},
0154 freq = S.f;
0155 xlbl_txt = 'Frequency [Hz]';
0156 ylbl1_txt = 'S(f) [m^2 s]';
0157 ylbl3_txt = 'Directional Spectrum';
0158 zlbl_txt = 'S(f,\theta) [m^2 s / rad]';
0159 funit = ' [Hz]';
0160 Sunit = ' [m^2 s ]';
0161 case {'k'},
0162 freq = S.k;
0163 xlbl_txt = 'Wave number [rad/m]';
0164 ylbl1_txt = 'S(k) [m^3/ rad]';
0165 funit = ' [rad/m]';
0166 Sunit = ' [m^3 / rad]';
0167 if isfield(S,'k2')
0168 ylbl4_txt='Wave Number Spectrum';
0169 end
0170 otherwise, error('Frequency type unknown')
0171 end
0172
0173 if isfield(S,'norm') & S.norm
0174 Sunit=[];funit=[];
0175 ylbl1_txt = 'Normalized Spectral density';
0176 ylbl3_txt = 'Normalized Directional Spectrum';
0177 ylbl4_txt = 'Normalized Wave Number Spectrum';
0178 if strcmp(ftype,'k')
0179 xlbl_txt = 'Normalized Wave number';
0180 else
0181 xlbl_txt = 'Normalized Frequency';
0182 end
0183 end
0184
0185 ylbl2_txt = 'Power spectrum (dB)';
0186
0187 phi =0;
0188 if isfield(S,'phi') & ~isempty(S.phi)
0189 phi = S.phi;
0190 end
0191
0192
0193 switch lower(S.type(end-2:end))
0194 case {'enc','req','k1d'}
0195 S.S = S.S(:);
0196 Fn = freq(end);
0197 indm = findpeaks(S.S,4);
0198 if isempty(indm)
0199
0200 [maxS,indm] = max(S.S);
0201 end
0202 maxS = max(S.S(indm));
0203 if isfield(S,'CI')& ~isempty(S.CI),
0204 maxS = maxS*S.CI(2);
0205 txtCI = [num2str(100*S.p), '% CI'];
0206 end
0207
0208 Fp = freq(indm);
0209
0210 if length(indm)==1
0211 txt = {['fp = ' num2str(Fp,2), funit];};
0212 else
0213 txt = {['fp1 = ' num2str(Fp(1),2), funit]};
0214 for j=2:length(indm)
0215 txt(end+1) = {['fp',num2str(j),' = ', num2str(Fp(j),2), funit]};
0216 end
0217 end
0218 if (plotflag == 3), subplot(2,1,1),end
0219 if (plotflag == 1) | (plotflag ==3),
0220 plot([Fp(:) Fp(:)]',[zeros(length(indm),1) S.S(indm)]',':')
0221 hold on
0222 if isfield(S,'CI'),
0223 plot(freq,S.S*S.CI(1), 'r:' )
0224 plot(freq,S.S*S.CI(2), 'r:' )
0225 end
0226 plot(freq,S.S,lintype) ,if ~ih, hold off,end
0227 if ih, a=axis; else a=zeros(1,4); end
0228 a1=Fn;
0229 if (Fp>0)
0230 a1=max(min(Fn,10*max(Fp)),a(2));
0231 end
0232 axis([0 a1 0 max(1.01*maxS,a(4))])
0233 title('Spectral density')
0234 xlabel(xlbl_txt)
0235 ylabel(ylbl1_txt )
0236 end
0237
0238 if (plotflag==3), subplot(2,1,2),end
0239
0240 if (plotflag == 2) | (plotflag ==3),
0241 ind=find(S.S>0);
0242
0243 plot([Fp(:),Fp(:)]',[repmat(min(10*log10(S.S(ind)/maxS)),length(Fp),1) 10*log10(S.S(indm)/maxS)]',':')
0244 hold on
0245 if isfield(S,'CI'),
0246 plot(freq(ind),10*log10(S.S(ind)*S.CI(1)/maxS), 'r:' )
0247 plot(freq(ind),10*log10(S.S(ind)*S.CI(2)/maxS), 'r:' )
0248 end
0249 plot(freq(ind),10*log10(S.S(ind)/maxS),lintype)
0250 if ~ih, hold off, end
0251 if ih, a=axis; else a=[0 0 0 0]; end
0252 axis([0 max(min(Fn,max(10*Fp)),a(2)) -20 max(1.01*10*log10(1),a(4))])
0253 title('Spectral density')
0254 xlabel(xlbl_txt)
0255 ylabel(ylbl2_txt )
0256 end
0257 if LegendOn
0258 if isfield(S,'CI'),
0259 legend(txt{:},txtCI,1)
0260 else
0261 legend(txt{:},1)
0262 end
0263 end
0264 case {'k2d'}
0265 if plotflag==1,
0266 [c, h] = contour(freq,S.k2,S.S,'b');
0267 z_level = clevels(c);
0268
0269
0270 textstart_x=0.05; textstart_y=0.94;
0271 cltext1(z_level,textstart_x,textstart_y);
0272 else
0273 [c,h] = contourf(freq,S.k2,S.S);
0274
0275 fcolorbar(c)
0276 end
0277 rotate(h,[0 0 1],-phi*180/pi)
0278
0279
0280
0281 xlabel(xlbl_txt)
0282 ylabel(xlbl_txt)
0283 title(ylbl4_txt)
0284
0285 km=max([-freq(1) freq(end) S.k2(1) -S.k2(end)]);
0286 axis([-km km -km km])
0287 hold on
0288 plot([0 0],[ -km km],':')
0289 plot([-km km],[0 0],':')
0290 axis('square')
0291
0292
0293
0294
0295 if ~ih, hold off,end
0296 case {'dir'}
0297 thmin = S.theta(1)-phi;thmax=S.theta(end)-phi;
0298 if plotflag==1
0299 if 0,
0300 h = polar([0 2*pi],[0 freq(end)]);
0301 delete(h);hold on
0302 [X,Y]=meshgrid(S.theta,freq);
0303 [X,Y]=pol2cart(X,Y);
0304 contour(X,Y,S.S',lintype)
0305 else
0306 if (abs(thmax-thmin)<3*pi),
0307 theta = S.theta;
0308 else
0309 theta = S.theta*pi/180;
0310 phi = phi*pi/180;
0311 end
0312 c = contours(theta,freq,S.S');
0313 if isempty(c)
0314 c = contours(theta,freq,S.S);
0315 end
0316 [z_level c] = clevels(c);
0317 h = polar(c(1,:),c(2,:),lintype);
0318 rotate(h,[0 0 1],-phi*180/pi)
0319 end
0320 title(ylbl3_txt)
0321
0322 textstart_x = -0.1; textstart_y=1.00;
0323 cltext1(z_level,textstart_x,textstart_y)
0324
0325 elseif (plotflag==2) | (plotflag==3),
0326
0327
0328 subplot(211)
0329
0330 if ih, hold on, end
0331
0332 Sf = spec2spec(S,'freq');
0333 wspecplot(Sf,1,lintype)
0334
0335 subplot(212)
0336
0337 Dtf = S.S;
0338 [Nt,Nf] = size(S.S);
0339 Sf = Sf.S(:).';
0340 ind = find(Sf);
0341
0342 if plotflag==3,
0343 Dtf(:,ind) = Dtf(:,ind)./Sf(ones(Nt,1),ind);
0344 end
0345 Dtheta = simpson(freq,Dtf,2);
0346 Dtheta = Dtheta/simpson(S.theta,Dtheta);
0347 [y,ind] = max(Dtheta);
0348 Wdir = S.theta(ind)-phi;
0349 txtwdir = ['\thetap=' num2pistr(Wdir,3)];
0350
0351 plot([1 1]*S.theta(ind)-phi,[0 Dtheta(ind)],':'), hold on
0352 if LegendOn
0353 lh=legend(txtwdir,0);
0354 end
0355 plot(S.theta-phi,Dtheta,lintype)
0356
0357 fixthetalabels(thmin,thmax,'x',2)
0358 ylabel('D(\theta)')
0359 title('Spreading function')
0360 if ~ih, hold off, end
0361
0362 elseif plotflag==4
0363 mesh(freq,S.theta-phi,S.S)
0364 xlabel(xlbl_txt);
0365 fixthetalabels(thmin,thmax,'y',3)
0366 zlabel(zlbl_txt)
0367 title(ylbl3_txt)
0368 elseif plotflag==5
0369
0370
0371 [X,Y]=meshgrid(S.theta-phi,freq);
0372 [X,Y]=pol2cart(X,Y);
0373 mesh(X,Y,S.S')
0374
0375 hold on, mesh(X,Y,zeros(size(S.S'))),hold off
0376 zlabel(zlbl_txt)
0377 title(ylbl3_txt)
0378 set(gca,'xticklabel','','yticklabel','')
0379 lighting phong
0380
0381
0382 elseif (plotflag==6) | (plotflag==7),
0383 theta = S.theta-phi;
0384 [c, h] = contour(freq,theta,S.S);
0385 fixthetalabels(thmin,thmax,'y',2)
0386 if plotflag==7,
0387 hold on
0388 [c,h] = contourf(freq,theta,S.S);
0389
0390 end
0391
0392 title(ylbl3_txt)
0393 xlabel(xlbl_txt);
0394 if 0,
0395 [z_level] = clevels(c);
0396
0397 textstart_x = 0.06; textstart_y=0.94;
0398 cltext1(z_level,textstart_x,textstart_y)
0399 else
0400 colormap('jet')
0401
0402 if plotflag==7,
0403 fcolorbar(c)
0404 else
0405
0406 hcb = colorbar;
0407 end
0408 grid on
0409 end
0410 else
0411 error('Unknown plot option')
0412 end
0413 otherwise, error('unknown spectral type')
0414 end
0415
0416 if ~ih, hold off, end
0417
0418
0419
0420
0421 set(findall(gcf,'type','text'),'buttondownfcn','edtext')
0422 set(gcf,'windowbuttondownfcn','edtext(''hide'')')
0423
0424 return
0425
0426 function xtxt = num2pistr(x,N)
0427
0428
0429
0430
0431
0432
0433
0434
0435
0436
0437
0438 if nargin<2|isempty(N)
0439 N=3;
0440 end
0441 [num den]=rat(x/pi);
0442 if (den<10)& (num<10)&(num~=0),
0443 if abs(den)==1, dtxt=[]; else, dtxt=['/' num2str(den)]; end
0444 if abs(num)==1,
0445 if num==-1, ntxt='-'; else, ntxt=[];end
0446 else
0447 ntxt=num2str(num);
0448 end
0449 xtxt= [ntxt '\pi' dtxt];
0450 else
0451 xtxt = num2str(x,N);
0452 end
0453 return
0454
0455 function fixthetalabels(thmin,thmax,xy,dim)
0456
0457
0458
0459
0460
0461
0462
0463
0464
0465
0466
0467 ind = [('x' == xy) ('y' == xy) ];
0468 yx = 'yx';
0469 yx = yx(ind);
0470 if nargin<4|isempty(dim),
0471 dim=2;
0472 end
0473
0474
0475
0476 if abs(thmax-thmin)<3*pi,
0477
0478 if xy=='x'
0479 if dim<3, axis([thmin,thmax 0 inf ]),else,axis([thmin,thmax 0 inf 0 inf]),end
0480 else
0481 if dim<3, axis([0 inf thmin,thmax ]),else,axis([0 inf thmin,thmax 0 inf]),end
0482 end
0483
0484 set(gca,[xy 'tick'],pi*(thmin/pi:0.25:thmax/pi));
0485 set(gca,[xy 'ticklabel'],[]);
0486 x = get(gca,[xy 'tick']);
0487 y = get(gca,[yx 'tick']);
0488 y1 = y(1);
0489 dy = y(2)-y1;
0490 yN = y(end)+dy;
0491 ylim = [y1 yN];
0492 dy1 = diff(ylim)/40;
0493
0494
0495 if xy=='x'
0496 for j=1:length(x)
0497 xtxt = num2pistr(x(j));
0498 figtext(x(j),y1-dy1,xtxt,'data','data','center','top');
0499 end
0500
0501 ax = [thmin thmax ylim];
0502 if dim<3,
0503 figtext(mean(x),y1-7*dy1,'Wave directions (rad)','data','data','center','top')
0504 else
0505 ax = [ax 0 inf];
0506 xlabel('Wave directions (rad)')
0507 end
0508 else
0509
0510 ax = [ylim thmin thmax];
0511
0512 if dim<3,
0513 for j=1:length(x)
0514 xtxt = num2pistr(x(j));
0515 figtext(y1-dy1/2,x(j),xtxt,'data','data','right');
0516 end
0517 set(gca,'DefaultTextRotation',90)
0518
0519 figtext(y1-3*dy1,mean(x),'Wave directions (rad)','data','data','center','bottom')
0520 set(gca,'DefaultTextRotation',0)
0521 else
0522 for j=1:length(x)
0523 xtxt = num2pistr(x(j));
0524 figtext(y1-3*dy1,x(j),xtxt,'data','data','right');
0525 end
0526 ax = [ax 0 inf];
0527 ylabel('Wave directions (rad)')
0528 end
0529 end
0530
0531
0532
0533
0534
0535 else
0536 set(gca,[xy 'tick'],thmin:45:thmax)
0537 if xy=='x'
0538 ax=[thmin thmax 0 inf];
0539 if dim>=3, ax=[ax 0 inf]; end
0540 xlabel('Wave directions (deg)')
0541 else
0542 ax=[0 inf thmin thmax ];
0543 if dim>=3, ax=[ax 0 inf]; end
0544 ylabel('Wave directions (deg)')
0545 end
0546 end
0547 axis(ax)
0548 return
0549
0550 function cltext1(z_level,textstart_x,textstart_y)
0551
0552
0553
0554
0555
0556
0557
0558
0559
0560
0561 if (nargin<4)|isempty(x_unit),
0562 x_unit='data';
0563 end
0564
0565
0566 delta_y = 1/33;
0567
0568
0569
0570
0571 h = figtext(textstart_x-0.05,textstart_y,' Level curves at:','norm');
0572 set(h,'FontWeight','Bold')
0573
0574
0575
0576
0577 textstart_y = textstart_y-delta_y;
0578 for ix=1:length(z_level)
0579 textstart_y = textstart_y-delta_y;
0580 figtext(textstart_x,textstart_y,num2str(z_level(ix),4),'norm');
0581 end
0582 return
0583