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
returns the frequency type of a Spectral density struct. | |
Get structure field contents. | |
Quick 1-D linear interpolation. | |
Linearly spaced vector. | |
Next higher power of 2. | |
Compare strings. |
Script to computer exercises 3 | |
Spectral simulation of a Gaussian sea, 2D (x,t) or 3D (x,y,t) | |
Computes covariance function and its derivatives, alternative version | |
Separates the linear component of the Spectrum | |
Simulates a Randomized 2nd order non-linear wave X(t) | |
Simulates a Gaussian process and its derivative from spectrum | |
Calculates joint density of minimum and following maximum |
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
Comments or corrections to the WAFO group