Home > wafo > multidim > private > bfsspec.m

bfsspec

PURPOSE ^

Estimate frequency spectrum for the surface elevation from the bfs timeseries

SYNOPSIS ^

SfBest = bfsSpec(Sf,Hw,pos,bfs);

DESCRIPTION ^

 BFSSPEC Estimate frequency spectrum for the surface elevation from the bfs timeseries 
  
   CALL:  SfBest = bfsSpec(Sf,Hw,pos,bfs);

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

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

Mathematical Statistics
Centre for Mathematical Sciences
Lund University with Lund Institute of Technology

Comments or corrections to the WAFO group


Generated on Thu 06-Oct-2005 02:21:16 for WAFO by m2html © 2003