Index: DPWavesProc/trunk/DPWavesProc/RELEASENOTES.txt =================================================================== --- adcp/trunk/adcp/RELEASENOTES.txt (revision 40) +++ DPWavesProc/trunk/DPWavesProc/RELEASENOTES.txt (revision 215) @@ -1,0 +1,1 @@ +This is the first BETA release of DPWavesProc. Index: DPWavesProc/trunk/DPWavesProc/adcp_matlab/radialtouvw.m =================================================================== --- adcp/trunk/adcp/adcp_matlab/radialtouvw.m (revision 168) +++ DPWavesProc/trunk/DPWavesProc/adcp_matlab/radialtouvw.m (revision 215) @@ -182,5 +182,5 @@ %transform the original x,y,z positions to the new positions accounting for %heading, pitch and roll... we also add adcpheight back in -newxyzpos=ones(3,4,5); +new_xyzpos=ones(3,4,5); new_xyzpos(1,:,:)=xyzpos(1,:,:)*a+xyzpos(2,:,:)*b+xyzpos(3,:,:)*c; new_xyzpos(2,:,:)=xyzpos(1,:,:)*d+xyzpos(2,:,:)*e+xyzpos(3,:,:)*f; Index: DPWavesProc/trunk/DPWavesProc/adcp_matlab/radialwaveplot.m =================================================================== --- adcp/trunk/adcp/adcp_matlab/radialwaveplot.m (revision 168) +++ DPWavesProc/trunk/DPWavesProc/adcp_matlab/radialwaveplot.m (revision 215) @@ -58,4 +58,17 @@ IMLMfreq=sum(real(radialI.S)')*dirres; wmonfreq=sum(wmon.S')*wmondirres; + +%Compute the coefficient for the upper and lower error bounds for the +%frequency +%spectrum assuming 95% confidence. Added on 9/17/08 +degF=radialE.degF; +chiUp=chi2inv(.975,degF); +chiLow=chi2inv(.025,degF); +coeffUp=degF/chiLow; +coeffLow=degF/chiUp; + +%calculate the conf limits throughout the frequency spectrum +EMEPfreqUP=EMEPfreq*coeffUp; +EMEPfreqLOW=EMEPfreq*coeffLow; %Find the maximum for the directional spectrum so we can set up the proper @@ -115,10 +128,12 @@ %plot the frequency energy spectrum subplot(1,2,2); -plot(freqs,EMEPfreq,'b'); +h1 = plot(freqs,EMEPfreq,'b'); hold on -plot(wmon.freqs,wmonfreq,'k'); +h2 = plot(freqs,EMEPfreqLOW,'b--'); +h3 = plot(freqs,EMEPfreqUP,'b--'); +h4 = plot(wmon.freqs,wmonfreq,'k'); axis(axis); -plot(freqs,IMLMfreq,'c'); -legend('EMEP radial vel','wavesmon','IMLM radial vel','location','best'); +%h5 = plot(freqs,IMLMfreq,'c'); +legend([h1,h2,h4],'EMEP radial vel','EMEP radial 95% conf limits', 'wavesmon','location','best'); title('directional wave spectrum integrated over direction'); xlabel('frequency in Hz'); @@ -135,4 +150,6 @@ % Calculate the Sig wave height EMEP_radial_Hsig=4*sqrt(m0); +%Use the function HsigConf.m to calculate the sigH confidence limits +EMEP_HsConf=HsigConf(radialE); % Calculate the peak period Tp @@ -151,4 +168,5 @@ disp(['EMEP radial vel']); disp(['SigH (meters): ' num2str(EMEP_radial_Hsig)]); +disp(['SigH 95% confidence limits: ' num2str(EMEP_HsConf)]); disp(['peak period (seconds): ' num2str(EMEP_radial_Tp)]); disp(['Dir of peak period: ' num2str(compangle(EMEP_radial_DTp, radialE.xaxisdir))]); @@ -157,4 +175,5 @@ E_radial_WI.hsig=EMEP_radial_Hsig; +E_radial_WI.hconf=EMEP_HsConf; E_radial_WI.tp=EMEP_radial_Tp; E_radial_WI.dtp=compangle(EMEP_radial_DTp, radialE.xaxisdir); @@ -169,4 +188,6 @@ % Calculate the Sig wave height IMLM_radial_Hsig=4*sqrt(m0); +%Use the function HsigConf.m to calculate the sigH confidence limits +IMLM_HsConf=HsigConf(radialI); % Calculate the peak period Tp @@ -185,4 +206,5 @@ disp(['IMLM radial vel']); disp(['SigH (meters): ' num2str(IMLM_radial_Hsig)]); +disp(['SigH 95% confidence limits: ' num2str(IMLM_HsConf)]); disp(['peak period (seconds): ' num2str(IMLM_radial_Tp)]); disp(['Dir of peak period: ' num2str(compangle(IMLM_radial_DTp, radialI.xaxisdir))]); @@ -191,4 +213,5 @@ I_radial_WI.hsig=IMLM_radial_Hsig; +I_radial_WI.hconf=IMLM_HsConf; I_radial_WI.tp=IMLM_radial_Tp; I_radial_WI.dtp=compangle(IMLM_radial_DTp, radialI.xaxisdir); Index: DPWavesProc/trunk/DPWavesProc/adcp_matlab/specmultiplot.m =================================================================== --- adcp/trunk/adcp/adcp_matlab/specmultiplot.m (revision 168) +++ DPWavesProc/trunk/DPWavesProc/adcp_matlab/specmultiplot.m (revision 215) @@ -35,4 +35,5 @@ Winfo.setup={'EMEP UVW','IMLM UVW','EMEP range','IMLM range','EMEP radial','IMLM radial','wavesmon'}'; Winfo.hsig=[]; +Winfo.hconf=[]; Winfo.peakP=[]; Winfo.dirP=[]; @@ -66,10 +67,10 @@ %run diwasp to generate this spectrum - [E_uvw_F01_D2, EPout]=dirspectest(ID,SM,EP,{'PLOTTYPE',0,'MESSAGE',0}); + [E_uvw_F01_D2, EPout]=dirspec(ID,SM,EP,{'PLOTTYPE',0,'MESSAGE',0}); %now for IMLM EP.method='IMLM'; EP.iter=3; - [I_uvw_F01_D2, EPout]=dirspectest(ID,SM,EP,{'PLOTTYPE',0,'MESSAGE',0}); + [I_uvw_F01_D2, EPout]=dirspec(ID,SM,EP,{'PLOTTYPE',0,'MESSAGE',0}); %set up data with ranges, freq and dir at default @@ -77,10 +78,10 @@ % run diwasp to generate this spectrum with EMEP - [E_range_F01_D2, EPout]=dirspectest(ID,SM,EP,{'PLOTTYPE',0,'MESSAGE',0}); + [E_range_F01_D2, EPout]=dirspec(ID,SM,EP,{'PLOTTYPE',0,'MESSAGE',0}); %now with IMLM EP.method='IMLM'; EP.iter=3; - [I_range_F01_D2, EPout]=dirspectest(ID,SM,EP,{'PLOTTYPE',0,'MESSAGE',0}); + [I_range_F01_D2, EPout]=dirspec(ID,SM,EP,{'PLOTTYPE',0,'MESSAGE',0}); %set up data with radial velocities, freq and dir at default @@ -88,10 +89,10 @@ % run diwasp to generate this spectrum with EMEP - [E_radial_F01_D2, EPout]=dirspectest(ID,SM,EP,{'PLOTTYPE',0,'MESSAGE',0}); + [E_radial_F01_D2, EPout]=dirspec(ID,SM,EP,{'PLOTTYPE',0,'MESSAGE',0}); %now with IMLM EP.method='IMLM'; EP.iter=3; - [I_radial_F01_D2, EPout]=dirspectest(ID,SM,EP,{'PLOTTYPE',0,'MESSAGE',0}); + [I_radial_F01_D2, EPout]=dirspec(ID,SM,EP,{'PLOTTYPE',0,'MESSAGE',0}); %Use waveplot to plot the polar, freq and dir plots of the spectra @@ -115,4 +116,5 @@ hsig=[vertcat(E_UVW_WI.hsig,I_UVW_WI.hsig,E_range_WI.hsig,I_range_WI.hsig,E_radial_WI.hsig,I_radial_WI.hsig,wmon_WI.hsig)]; + hconf=[vertcat(E_UVW_WI.hconf,I_UVW_WI.hconf,E_range_WI.hconf,I_range_WI.hconf,E_radial_WI.hconf,I_radial_WI.hconf,[NaN NaN])]; peakP=[vertcat(E_UVW_WI.tp,I_UVW_WI.tp,E_range_WI.tp,I_range_WI.tp,E_radial_WI.tp,I_radial_WI.tp,wmon_WI.tp)]; dirP=[vertcat(E_UVW_WI.dtp,I_UVW_WI.dtp,E_range_WI.dtp,I_range_WI.dtp,E_radial_WI.dtp,I_radial_WI.dtp,wmon_WI.dtp)]; @@ -121,4 +123,5 @@ Winfo.hsig=[horzcat(Winfo.hsig,hsig)]; + Winfo.hconf=[horzcat(Winfo.hconf,hconf)]; Winfo.peakP=[horzcat(Winfo.peakP,peakP)]; Winfo.dirP=[horzcat(Winfo.dirP,dirP)]; Index: DPWavesProc/trunk/DPWavesProc/adcp_matlab/waveplot.m =================================================================== --- adcp/trunk/adcp/adcp_matlab/waveplot.m (revision 168) +++ DPWavesProc/trunk/DPWavesProc/adcp_matlab/waveplot.m (revision 215) @@ -70,4 +70,12 @@ IMLMrangefreq=sum(real(rangeI.S)')*dirres; wmonfreq=sum(wmon.S')*wmondirres; + +%Compute the coefficient for the upper and lower error bounds for the power +%spectrum assuming 95% confidence. Added on 9/17/08 +degF=EMEP.degF; +chiUp=chi2inv(.975,degF); +chiLow=chi2inv(.025,degF); +coeffUp=degF/chiLow; +coeffLow=degF/chiUp; %Find the maximum for the directional spectrum so we can set up the proper @@ -148,11 +156,8 @@ %For EMEP uvw - -%calculate the 0,1,2 moments -m0=sum(EMEPfreq*freqres); -m1=sum(freqs.*EMEPfreq*freqres); -m2=sum((freqs.^2).*EMEPfreq*freqres); -% Calculate the Sig wave height -EMEP_UVW_Hsig=4*sqrt(m0); +% Calculate the Sig wave height +EMEP_UVW_Hsig=Hsig(EMEP); +%Use the function HsigConf.m to calculate the sigH confidence limits +EMEP_HsConf=HsigConf(EMEP); % Calculate the peak period Tp @@ -171,4 +176,5 @@ disp(['EMEP uvw']); disp(['SigH (meters): ' num2str(EMEP_UVW_Hsig)]); +disp(['SigH 95% confidence limits: ' num2str(EMEP_HsConf)]); disp(['peak period (seconds): ' num2str(EMEP_UVW_Tp)]); disp(['Dir of peak period: ' num2str(compangle(EMEP_UVW_DTp, EMEP.xaxisdir))]); @@ -177,4 +183,5 @@ E_UVW_WI.hsig=EMEP_UVW_Hsig; +E_UVW_WI.hconf=EMEP_HsConf; E_UVW_WI.tp=EMEP_UVW_Tp; E_UVW_WI.dtp=compangle(EMEP_UVW_DTp, EMEP.xaxisdir); @@ -184,10 +191,8 @@ % For IMLM uvw -%calculate the 0,1,2 moments -m0=sum(IMLMfreq*freqres); -m1=sum(freqs.*IMLMfreq*freqres); -m2=sum((freqs.^2).*IMLMfreq*freqres); -% Calculate the Sig wave height -IMLM_UVW_Hsig=4*sqrt(m0); +% Calculate the Sig wave height +IMLM_UVW_Hsig=Hsig(IMLM); +%Use the function HsigConf.m to calculate the sigH confidence limits +IMLM_HsConf=HsigConf(IMLM); % Calculate the peak period Tp @@ -206,4 +211,5 @@ disp(['IMLM uvw']); disp(['SigH (meters): ' num2str(IMLM_UVW_Hsig)]); +disp(['SigH 95% confidence limits: ' num2str(IMLM_HsConf)]); disp(['peak period (seconds): ' num2str(IMLM_UVW_Tp)]); disp(['Dir of peak period: ' num2str(compangle(IMLM_UVW_DTp, IMLM.xaxisdir))]); @@ -212,4 +218,5 @@ I_UVW_WI.hsig=IMLM_UVW_Hsig; +I_UVW_WI.hconf=IMLM_HsConf; I_UVW_WI.tp=IMLM_UVW_Tp; I_UVW_WI.dtp=compangle(IMLM_UVW_DTp, IMLM.xaxisdir); @@ -218,10 +225,8 @@ % for Range of EMEP -%calculate the 0,1,2 moments -m0=sum(EMEPrangefreq*freqres); -m1=sum(freqs.*EMEPrangefreq*freqres); -m2=sum((freqs.^2).*EMEPrangefreq*freqres); -% Calculate the Sig wave height -EMEP_range_Hsig=4*sqrt(m0); +% Calculate the Sig wave height +EMEP_range_Hsig=Hsig(rangeE); +%Use the function HsigConf.m to calculate the sigH confidence limits +EMEP_range_HsConf=HsigConf(rangeE); % Calculate the peak period Tp @@ -240,4 +245,5 @@ disp(['EMEP range']); disp(['SigH (meters): ' num2str(EMEP_range_Hsig)]); +disp(['SigH 95% confidence limits: ' num2str(EMEP_range_HsConf)]); disp(['peak period (seconds): ' num2str(EMEP_range_Tp)]); disp(['Dir of peak period: ' num2str(compangle(EMEP_range_DTp, rangeE.xaxisdir))]); @@ -246,4 +252,5 @@ E_range_WI.hsig=EMEP_range_Hsig; +E_range_WI.hconf=EMEP_range_HsConf; E_range_WI.tp=EMEP_range_Tp; E_range_WI.dtp=compangle(EMEP_range_DTp, rangeE.xaxisdir); @@ -253,10 +260,8 @@ %for Range of IMLM -%calculate the 0,1,2 moments -m0=sum(IMLMrangefreq*freqres); -m1=sum(freqs.*IMLMrangefreq*freqres); -m2=sum((freqs.^2).*IMLMrangefreq*freqres); -% Calculate the Sig wave height -IMLM_range_Hsig=4*sqrt(m0); +% Calculate the Sig wave height +IMLM_range_Hsig=Hsig(rangeI); +%Use the function HsigConf.m to calculate the sigH confidence limits +IMLM_range_HsConf=HsigConf(rangeI); % Calculate the peak period Tp @@ -275,4 +280,5 @@ disp(['IMLM range']); disp(['SigH (meters): ' num2str(IMLM_range_Hsig)]); +disp(['SigH 95% confidence limits: ' num2str(IMLM_range_HsConf)]); disp(['peak period (seconds): ' num2str(IMLM_range_Tp)]); disp(['Dir of peak period: ' num2str(compangle(IMLM_range_DTp, rangeI.xaxisdir))]); @@ -281,4 +287,5 @@ I_range_WI.hsig=IMLM_range_Hsig; +I_range_WI.hconf=IMLM_range_HsConf; I_range_WI.tp=IMLM_range_Tp; I_range_WI.dtp=compangle(IMLM_range_DTp, rangeI.xaxisdir); Index: DPWavesProc/trunk/DPWavesProc/diwasp_1_1GD/dirspec.m =================================================================== --- adcp/trunk/adcp/diwasp_1_1GD/dirspec.m (revision 168) +++ DPWavesProc/trunk/DPWavesProc/diwasp_1_1GD/dirspec.m (revision 215) @@ -1,3 +1,3 @@ -function [SM,EP]=dirspec(ID,SM,EP,varargin); +function [SM,EP]=dirspec(ID,SM,EP,varargin) %DIWASP V1.1 function with GD edit @@ -84,11 +84,29 @@ end -% NOTE: changed from csd to cpsd 6/15/07 -%calculate the cross-power spectra +%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 - [xpstmp,Ftmp]=cpsd(data(:,m),data(:,n),[],[],nfft,ID.fs); + % 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); - end +end end @@ -119,12 +137,29 @@ end +%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 compare=strcmp('radial',ID.datatypes(1)); - if compare == 1 Ss=mean(Ss); -end - + %since this is using radial velocities the number of data series used + %the total number of sensors or szd + dataUsed=szd; +end + + + +%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=dataUsed*2*numWindows*windowConstant; +SM.degF=degF; + + ffs=sum(F<=max(SM.freqs)); Index: DPWavesProc/trunk/DPWavesProc/diwasp_1_1GD/infospec.m =================================================================== --- adcp/trunk/adcp/diwasp_1_1GD/infospec.m (revision 168) +++ DPWavesProc/trunk/DPWavesProc/diwasp_1_1GD/infospec.m (revision 215) @@ -1,3 +1,6 @@ -function [H,Tp,DTp,Dp]=infospec(SM) +function [H,HsConf,Tp,DTp,Dp]=infospec(SM) + +%Edited on 9/17/08 to include confidence limits for the sig wave height + %DIWASP V1.1 function @@ -8,4 +11,5 @@ %Outputs: %Hsig Signficant wave height (Hmo) +%HsConf 95% conf limits for Hsig %Tp Peak period %DTp Direction of spectral peak @@ -14,4 +18,5 @@ %Inputs: %SM A spectral matrix structure containing the file data +% % %Hsig is the significant wave height. Tp is the peak frequency, the highest point in the one dimensional spectrum. @@ -25,5 +30,6 @@ SM=check_data(SM,2);if isempty(SM) return;end; -H=HSig(SM); +H=Hsig(SM); +HsConf=HsigConf(SM); S=sum(real(SM.S),2); % S is freq spectrum @@ -38,4 +44,5 @@ disp(['Infospec::']); disp(['Significant wave height (Hmo): ' num2str(H)]); +disp(['Sig wave height 95% confidence limits: ' num2str(HsConf)]); disp(['Peak period: ' num2str(Tp)]); disp(['Direction of peak period: ' num2str(DTp) ' axis angle / ' num2str(compangle(DTp,SM.xaxisdir)) ' compass bearing']); Index: DPWavesProc/trunk/DPWavesProc/diwasp_1_1GD/private/Hsig.m =================================================================== --- adcp/trunk/adcp/diwasp_1_1GD/private/Hsig.m (revision 168) +++ DPWavesProc/trunk/DPWavesProc/diwasp_1_1GD/private/Hsig.m (revision 215) @@ -1,3 +1,3 @@ -function Hs=Hmo(varargin) +function Hs=Hsig(varargin) %DIWASP V1.1 function to calculate significant wave height % Index: DPWavesProc/trunk/DPWavesProc/nortek/nortek_Winfo_plot.m =================================================================== --- adcp/trunk/adcp/nortek/nortek_Winfo_plot.m (revision 168) +++ DPWavesProc/trunk/DPWavesProc/nortek/nortek_Winfo_plot.m (revision 215) @@ -10,9 +10,13 @@ h1=plot(Winfo.time,Winfo.hsig(1,:),'b.-'); hold on +h1a=plot(Winfo.time,Winfo.hconf(1,1:2:end),'b--'); +plot(Winfo.time,Winfo.hconf(1,2:2:end),'b--'); h2=plot(Winfo.time,Winfo.hsig(2,:),'r.-'); +h2a=plot(Winfo.time,Winfo.hconf(2,1:2:end),'r--'); +plot(Winfo.time,Winfo.hconf(2,2:2:end),'r--'); h3=plot(Winfo.time,Winfo.hsig(3,:),'k.-'); datetick('x',15,'keepticks'); -legend([h1, h2, h3],'EMEP','IMLM','nortek','location','best'); +legend([h1,h1a,h2,h2a,h3],'EMEP','EMEP 95% confidence','IMLM','IMLM 95% confidence','nortek','location','best'); title(['Significant Wave Height starting on ',datestr(Winfo.time(:,1),0)]); xlabel('time'); Index: DPWavesProc/trunk/DPWavesProc/nortek/nortek_specmultiplot.m =================================================================== --- adcp/trunk/adcp/nortek/nortek_specmultiplot.m (revision 168) +++ DPWavesProc/trunk/DPWavesProc/nortek/nortek_specmultiplot.m (revision 215) @@ -11,5 +11,5 @@ %peak period, direction of peak period and dominant dir. It also generates %polar plots and 2d freq and directional plots of a variety of different -%estimation methods, resolutions and the wavesmon output. +%estimation methods, resolutions and the nortek quickwave output. @@ -38,4 +38,5 @@ Winfo.setup={'EMEP','IMLM' 'nortek'}'; Winfo.hsig=[]; +Winfo.hconf=[]; Winfo.peakP=[]; Winfo.dirP=[]; @@ -68,10 +69,10 @@ %run diwasp to generate this spectrum - [E_F01_D2, EPout]=dirspectest(ID,SM,EP,{'PLOTTYPE',0,'MESSAGE',0}); + [E_F01_D2, EPout]=dirspec(ID,SM,EP,{'PLOTTYPE',0,'MESSAGE',0}); %now for IMLM EP.method='IMLM'; EP.iter=3; - [I_F01_D2, EPout]=dirspectest(ID,SM,EP,{'PLOTTYPE',0,'MESSAGE',0}); + [I_F01_D2, EPout]=dirspec(ID,SM,EP,{'PLOTTYPE',0,'MESSAGE',0}); @@ -82,9 +83,11 @@ saveas(fig1,['polar1_' dateinfo '.fig']); saveas(fig2,['dirfreq1_' dateinfo '.fig']); + + close all; - %Fill up the Winfo data structure hsig=[vertcat(E_WI.hsig,I_WI.hsig,nortek_WI.hsig)]; + hconf=[vertcat(E_WI.hconf,I_WI.hconf,nortek_WI.hconf)]; peakP=[vertcat(E_WI.tp,I_WI.tp,nortek_WI.tp)]; dirP=[vertcat(E_WI.dtp,I_WI.dtp,nortek_WI.dtp)]; @@ -92,4 +95,5 @@ Winfo.hsig=[horzcat(Winfo.hsig,hsig)]; + Winfo.hconf=[horzcat(Winfo.hconf,hconf)]; Winfo.peakP=[horzcat(Winfo.peakP,peakP)]; Winfo.dirP=[horzcat(Winfo.dirP,dirP)]; Index: DPWavesProc/trunk/DPWavesProc/nortek/nortek_waveplot.m =================================================================== --- adcp/trunk/adcp/nortek/nortek_waveplot.m (revision 168) +++ DPWavesProc/trunk/DPWavesProc/nortek/nortek_waveplot.m (revision 215) @@ -54,4 +54,19 @@ IMLMfreq=sum(real(IMLM.S)')*dirres; nortekfreq=sum(nortek.S')*nortekdirres; + +%Compute the coefficient for the upper and lower error bounds for the +%frequency +%spectrum assuming 95% confidence. Added on 9/17/08 +degF=2*EMEP.degF; % 2* since for the Nortek the range beam sampled at 4hz (4096) is decimated or averaged to 2hz (2048) +chiUp=chi2inv(.975,degF); +chiLow=chi2inv(.025,degF); +coeffUp=degF/chiLow; +coeffLow=degF/chiUp; + +%calculate the conf limits throughout the frequency spectrum +EMEPfreqUP=EMEPfreq*coeffUp; +EMEPfreqLOW=EMEPfreq*coeffLow; +IMLMfreqUP=IMLMfreq*coeffUp; +IMLMfreqLOW=IMLMfreq*coeffLow; %Find the maximum for the directional spectrum so we can set up the proper @@ -110,9 +125,13 @@ %plot the frequency energy spectrum subplot(1,2,2); -plot(freqs,EMEPfreq,'b'); +h1=plot(freqs,EMEPfreq,'b'); hold on -plot(freqs,IMLMfreq,'r'); -plot(nortekfreqs,nortekfreq,'k'); -legend('EMEP','IMLM','nortek quickwave','location','best'); +h2=plot(freqs,EMEPfreqUP,'b--'); +plot(freqs,EMEPfreqLOW,'b--'); +h3=plot(freqs,IMLMfreq,'r'); +h4=plot(freqs,IMLMfreqUP,'r--'); +plot(freqs,IMLMfreqLOW,'r--'); +h5=plot(nortekfreqs,nortekfreq,'k'); +legend('EMEP','EMEP 95% confidence','IMLM','IMLM 95% confidence','nortek quickwave','location','best'); title('directional wave spectrum integrated over direction'); xlabel('frequency in Hz'); @@ -129,4 +148,6 @@ % Calculate the Sig wave height EMEP_Hsig=4*sqrt(m0); +%Use the function HsigConf.m to calculate the sigH confidence limits +EMEP_HsConf=nortek_HsigConf(EMEP); % Calculate the peak period Tp @@ -145,4 +166,5 @@ disp(['EMEP']); disp(['SigH (meters): ' num2str(EMEP_Hsig)]); +disp(['SigH 95% conf. (meters): ' num2str(EMEP_HsConf)]); disp(['peak period (seconds): ' num2str(EMEP_Tp)]); disp(['Dir of peak period: ' num2str(compangle(EMEP_DTp, EMEP.xaxisdir))]); @@ -151,4 +173,5 @@ E_WI.hsig=EMEP_Hsig; +E_WI.hconf=EMEP_HsConf; E_WI.tp=EMEP_Tp; E_WI.dtp=compangle(EMEP_DTp, EMEP.xaxisdir); @@ -164,4 +187,6 @@ % Calculate the Sig wave height IMLM_Hsig=4*sqrt(m0); +%Use the function HsigConf.m to calculate the sigH confidence limits +IMLM_HsConf=nortek_HsigConf(IMLM); % Calculate the peak period Tp @@ -180,4 +205,5 @@ disp(['IMLM']); disp(['SigH (meters): ' num2str(IMLM_Hsig)]); +disp(['SigH 95% conf. (meters): ' num2str(IMLM_HsConf)]); disp(['peak period (seconds): ' num2str(IMLM_Tp)]); disp(['Dir of peak period: ' num2str(compangle(IMLM_DTp, IMLM.xaxisdir))]); @@ -186,4 +212,5 @@ I_WI.hsig=IMLM_Hsig; +I_WI.hconf=IMLM_HsConf; I_WI.tp=IMLM_Tp; I_WI.dtp=compangle(IMLM_DTp, IMLM.xaxisdir); @@ -222,4 +249,5 @@ nortek_WI.hsig=nortek_Hsig; +nortek_WI.hconf=[NaN NaN]; nortek_WI.tp=nortek_Tp; nortek_WI.dtp=compangle(nortek_DTp, nortek.xaxisdir); Index: DPWavesProc/trunk/DPWavesProc/radialtouvw1.m =================================================================== --- adcp/trunk/adcp/radialtouvw1.m (revision 168) +++ DPWavesProc/trunk/DPWavesProc/radialtouvw1.m (revision 215) @@ -1,3 +1,3 @@ -function [ID,SM,EP,orbitnew]=radialtouvw(pressure,range,orbit,sysinfo,data_type) +function [ID,SM,EP,orbitnew]=radialtouvw1(pressure,range,orbit,sysinfo,data_type)