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:
 error Display message and abort function. mean Average or mean value. unique Set unique.
This function is called by:
 dat2dspec Estimates the directional wave spectrum from timeseries

## 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