BFSSPEC Estimate frequency spectrum for the surface elevation from the bfs timeseries CALL: SfBest = bfsSpec(Sf,Hw,pos,bfs);
Estimates the directional wave spectrum from timeseries |
001 function SfBest = bfsSpec(Sf,Hw,pos,bfs); 002 %BFSSPEC Estimate frequency spectrum for the surface elevation from the bfs timeseries 003 % 004 % CALL: SfBest = bfsSpec(Sf,Hw,pos,bfs); 005 % 006 % 007 def = unique(pos(bfs,4)).'; 008 009 k0 = find(def==5|def==7|def==11|def==14|def==17); 010 if any(k0), def(k0) = def(k0)-1;end 011 012 k0 = find(def==8); 013 if any(k0), def(k0) = def(k0)-2;end 014 015 def = unique(def); 016 SfBest = Sf; 017 018 for sensorType1 = def 019 if (sensorType1==6) % surface curvatures : n_xx, n_yy, n_xy 020 sensorType2 = 7; sensorType3 = 8; 021 022 %Find the indices to sensorTypes 1,2 and 3: 023 ix1 = find(pos(:,4)==sensorType1); Nx1 = length(Nx1); 024 ix2 = find(pos(:,4)==sensorType2); Nx2 = length(Nx2); 025 ix3 = find(pos(:,4)==sensorType3); Nx3 = length(Nx3); 026 027 Nx = min([Nx1,Nx2,Nx3]); % need at least one pair of observations 028 if (Nx>0) 029 k0 = find(Hw(ix1(1),:)==0); 030 if any(k0),SfBest([ix1;ix2;ix3],k0)=0 ;end 031 032 k = find(Hw(ix1(1),:)~=0); 033 Sf0 = mean(Sf(ix1(1:Nx),k)+Sf(ix2(1:Nx),k)+2*Sf(ix3(1:Nx),k),1)./(Hw(ix1(1),k).^2); 034 SfBest([ix1;ix2;ix3],k) = Sf0(ones(Nx1+Nx2+Nx3,1),k); 035 else 036 error('Unable to estimate the sea surface spectrum from the surface curvature!') 037 end 038 039 elseif any(sensorType1==[2:3 9 12 15 18]) 040 041 % Find indices to sensorType1 042 ix1 = find(pos(:,4)==sensorType1); 043 Nx = length(ix1); 044 if Nx>0 045 k0 = find(Hw(ix1(1),:)==0); 046 if any(k0),SfBest(ix1,k0)=0 ;end 047 048 k = find(Hw(ix1(1),:)~=0); 049 Sf0 = mean(Sf(ix1,k),1)./(Hw(ix1(1),k).^2); 050 SfBest(ix1,k) = Sf0(ones(Nx,1),:); 051 else 052 error('Unable to estimate the sea surface spectrum') 053 end 054 elseif any(sensorType1==[4 10 14 16]) 055 % surface slopes,water particle velocity, water particle acceleration 056 % or water particle displacement 057 sensorType2 = sensorType1+1; 058 059 % Find indices to sensorType1 and sensorType2 060 ix1 = find(pos(:,4)==sensorType1); Nx1 = length(ix1); 061 ix2 = find(pos(:,4)==sensorType2); Nx2 = length(ix2); 062 Nx = min(Nx1,Nx2); % need at least one pair of observations 063 if (Nx>0) 064 k0 = find(Hw(ix1(1),:)==0); 065 if any(k0),SfBest([ix1;ix2],k0)=0 ;end 066 067 k = find(Hw(ix1(1),:)~=0); 068 Sf0 = mean(Sf(ix1(1:Nx),k)+Sf(ix2(1:Nx),k),1)./(Hw(ix1(1),k).^2); 069 SfBest([ix1;ix2],k) = Sf0(ones(Nx1+Nx2,1),:); 070 else %if (any() 071 error('Unable to estimate the sea surface spectrum') 072 end 073 end 074 end 075 076 %Keep only the best frequency spectra 077 SfBest = SfBest(bfs,:); 078 mmx = max(SfBest(:))*1e-5; % Minimum value which is not considered as noise! 079 k = find(SfBest < mmx); % This is most probably noise 080 if 0, % Geometric mean 081 if any(k) 082 SfBest(k) = abs(SfBest(k)+mmx*1e-10); 083 end 084 if length(bfs) > 1, 085 SfBest = exp(mean(log(SfBest))); % geometric mean 086 end 087 else % Ordinary mean 088 if any(k), % Set the noise to zero. 089 SfBest(k) = 0; 090 end 091 if length(bfs) > 1, 092 SfBest = mean(SfBest); % ordinary mean 093 end 094 end 095 return
Comments or corrections to the WAFO group