Home > wafo > spec > specinterp.m

specinterp

PURPOSE ^

Interpolation and zero-padding of spectrum

SYNOPSIS ^

Snew = specinterp(S,dt)

DESCRIPTION ^

 SPECINTERP Interpolation and zero-padding of spectrum
             to change Nyquist freq.
 
  CALL:  Snew = specinterp(S,dt)
 
     Snew, S = spectrum structs (any type)
          dt = wanted sampling interval (default as given by S, see spec2dt)
               unit: [s] if frequency-spectrum, [m] if wave number spectrum  
               
  To be used before simulation (e.g. spec2sdat) or evaluation of covariance
  function (spec2cov) to directly get wanted sampling interval.
  The input spectrum is interpolated and padded with zeros to reach
  the right max-frequency, w(end)=pi/dt, f(end)=1/2/dt, or k(end)=pi/dt.
  The objective is that output frequency grid should be at least as dense as
  the input grid, but the number of points in the spectrum is maximized to
  2^13+1. 
  NB! Also zero-padding down to zero freq, if S does not start there. 
      If empty input dt, this is the only effect.
   
  See also  spec2cov, spec2sdat, covinterp, spec2dt

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

001 function Snew = specinterp(S,dt)
002 %SPECINTERP Interpolation and zero-padding of spectrum
003 %            to change Nyquist freq.
004 %
005 % CALL:  Snew = specinterp(S,dt)
006 %
007 %    Snew, S = spectrum structs (any type)
008 %         dt = wanted sampling interval (default as given by S, see spec2dt)
009 %              unit: [s] if frequency-spectrum, [m] if wave number spectrum  
010 %              
011 % To be used before simulation (e.g. spec2sdat) or evaluation of covariance
012 % function (spec2cov) to directly get wanted sampling interval.
013 % The input spectrum is interpolated and padded with zeros to reach
014 % the right max-frequency, w(end)=pi/dt, f(end)=1/2/dt, or k(end)=pi/dt.
015 % The objective is that output frequency grid should be at least as dense as
016 % the input grid, but the number of points in the spectrum is maximized to
017 % 2^13+1. 
018 % NB! Also zero-padding down to zero freq, if S does not start there. 
019 %     If empty input dt, this is the only effect.
020 %  
021 % See also  spec2cov, spec2sdat, covinterp, spec2dt
022 
023 % Tested on:
024 % History:
025 % revised pab 12.10.2001
026 % -fixed a bug created 11.10.2001: ftype='k' now works OK
027 % revised pab 11.10.2001
028 % - added call to freqtype.m
029 % - fixed a bug: ftype=='f' now works correctly
030 % revised es 080600, revision of last revision, output matrix now OK
031 %  revised by es 22.05.00, output vectors columns, not rows
032 %  by es 13.01.2000, original ideas by sylvie, ir
033   
034 Snew  = S;  
035 ftype = freqtype(S);
036 w     = getfield(S,ftype);
037 n     = length(w);
038 
039 if strcmp(ftype,'f') %ftype==f
040   dT=1/(2*w(n));% sampling interval=1/Fs
041 else % ftype == w og ftype == k
042   dT=pi/w(n);
043 end
044 if nargin<2|isempty(dt),  dt=dT; end
045 
046 % Find how many points that is needed
047 nfft   = 2^nextpow2(n-1);
048 dttest = dT*(n-1)/nfft;
049 while (dttest>dt) & (nfft<2^13)
050  nfft=nfft*2;
051  dttest=dT*(n-1)/nfft;
052 end;
053 nfft=nfft+1;
054 
055 if strcmp(ftype,'w')
056   Snew.w=linspace(0,pi/dt,nfft)';
057 elseif strcmp(ftype,'f') %freqtype==f
058   Snew.f=linspace(0,1/2/dt,nfft)';
059 else
060   Snew.k=linspace(0,pi/dt,nfft)';
061 end
062 
063 if size(S.S,2)<2
064   S.S=S.S'; % if vector, make it a row to match matrix (np x nf)
065 end
066 w=w(:);
067 dw=w(end)-w(end-1);
068 
069 fmax = max(getfield(Snew,ftype));
070 if fmax>w(end)
071   % add a zero just above old max-freq, and a zero at new max-freq
072   % to get correct interpolation there
073   w=[w;w(end)+dw; fmax ];
074   S.S=[S.S zeros(size(S.S,1),2)];
075 end
076 if w(1)>0
077   % add a zero at freq 0, and, if there is space, a zero just below min-freq
078   if w(1)>dw
079     w=[w(1)-dw;w];
080     S.S=[zeros(size(S.S,1),1) S.S];
081   end
082   size(w)
083   w=[0;w];
084   S.S=[zeros(size(S.S,1),1) S.S];
085 end  
086 
087 Snew.S=interp1q(w,S.S',getfield(Snew,ftype))';
088 
089 if size(Snew.S,1)<2
090   Snew.S=Snew.S.'; % if vector, make it a column
091 end
092 
093 
094 
095 
096

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