COV2SPEC Computes spectral density given the auto covariance function CALL: S = cov2spec(R,rate,nugget,trunc,fast); R = a covariance structure S = a spectral density structure Optional parameters: rate = 1,2,4,8...2^r, interpolation rate for f (default = 1, no interpolation) nugget = add a nugget effect to ensure that round off errors do not result in negative spectral estimates (default 0) Good choice might be 10^-12. trunc = truncates all spectral values where S/max(S) < trunc 0 <= trunc <1 This is to ensure that high frequency noise is not added to the spectrum. (default 1e-5) fast = 1 if zero-pad to obtain power of 2 length R.R (default) 0 if no zero-padding of R.R, slower but more accurate. NB! This routine requires that the covariance is evenly spaced starting from zero lag. Currently only capable of 1D matrices. Example: R = createcov; L = 129; R.t = linspace(0,75,L).'; R.R = zeros(L,1); win = parzen(40); R.R(1:20) = win(21:40); S = cov2spec(R); See also spec2cov, datastructures
Spectrum structure constructor | |
Display message and abort function. | |
Discrete Fourier transform. | |
Get structure field names. | |
Get structure field contents. | |
Resample data at a higher rate using lowpass interpolation. | |
True for structures. | |
Linearly spaced vector. | |
Convert string to lowercase. | |
Next higher power of 2. | |
Set structure field contents. | |
Compare strings. | |
Compare strings ignoring case. |
001 function S = cov2spec(R,rate,nugget,trunc,fast); 002 %COV2SPEC Computes spectral density given the auto covariance function 003 % 004 % CALL: S = cov2spec(R,rate,nugget,trunc,fast); 005 % 006 % R = a covariance structure 007 % 008 % S = a spectral density structure 009 % 010 % Optional parameters: 011 % 012 % rate = 1,2,4,8...2^r, interpolation rate for f 013 % (default = 1, no interpolation) 014 % 015 % nugget = add a nugget effect to ensure that round off errors 016 % do not result in negative spectral estimates 017 % (default 0) Good choice might be 10^-12. 018 % 019 % trunc = truncates all spectral values where S/max(S) < trunc 020 % 0 <= trunc <1 This is to ensure that high frequency 021 % noise is not added to the spectrum. (default 1e-5) 022 % 023 % fast = 1 if zero-pad to obtain power of 2 length R.R (default) 024 % 0 if no zero-padding of R.R, slower but more accurate. 025 % 026 % NB! This routine requires that the covariance is evenly spaced 027 % starting from zero lag. Currently only capable of 1D matrices. 028 % 029 % Example: 030 % R = createcov; 031 % L = 129; 032 % R.t = linspace(0,75,L).'; 033 % R.R = zeros(L,1); 034 % win = parzen(40); 035 % R.R(1:20) = win(21:40); 036 % S = cov2spec(R); 037 % 038 % See also spec2cov, datastructures 039 040 % tested on: matlab 5.3 041 % history: 042 % Revised pab 043 % -Increased the precision in calculations 044 % - deleted L from input, added fast instead 045 % revised jr 11.07.2000 046 % - line 158. Changed n to nf . 047 % revised pab 13.03.2000 048 % - commented out warning messages 049 % revised pab 16.01.2000 050 % - fixed a bug: truncate negative values of the spectral density to 051 % zero to make sure these spurious points are not added to the 052 % spectrum 053 % - added trunc: this is a fix: truncating the values of the spectral density to 054 % zero if S/max(S) < trunc. This is to ensure that high frequency 055 % noise is not added to the spectrum. 056 % revised by pab 23.09.1999 057 058 % dT = time-step between data points.(default = T(2)-T1)). 059 060 % add a nugget effect to ensure that round off errors 061 % do not result in negative spectral estimates 062 063 if ~isstruct(R) 064 error('Incorrect input Covariance, see help datastructures') 065 end 066 if ndims(R.R)>2|(min(size(R.R))>1), 067 error('This function is only capable of 1D covariances') 068 end 069 070 if nargin<3|isempty(nugget) 071 nugget = 0;%10^-12; 072 end 073 074 if nargin<4|isempty(trunc) 075 trunc = 1e-5; 076 else 077 trunc = min(abs(trunc),1); % make sure it 078 end 079 if nargin<5 | isempty(fast) 080 fast = 1; 081 end 082 083 comment=0; 084 085 086 087 n = length(R.R); 088 names = fieldnames(R); 089 ind = find(strcmpi(names,'x')+strcmp(names,'y')+strcmp(names,'t')); 090 %options are 'x' and 'y' and 't' 091 092 lagtype = lower(names{ind}); 093 if strcmpi(lagtype,'t') 094 spectype = 'freq'; 095 ftype = 'w'; %options f=S(f), w=S(w) 096 else 097 spectype = 'k1d'; 098 ftype = 'k'; %options f=S(f), w=S(w) 099 end 100 ti = getfield(R,lagtype); 101 dT = ti(2)-ti(1); % ti 102 103 104 105 if nargin<2|isempty(rate) 106 rate = 1;%interpolation rate 107 else 108 rate = 2^nextpow2(rate);%make sure rate is a power of 2 109 end 110 111 % add a nugget effect to ensure that round off errors 112 % do not result in negative spectral estimates 113 ACF = R.R(:); 114 ACF(1) = R.R(1) +nugget; 115 % embedding a circulant vector and Fourier transform 116 if fast 117 nfft = 2^nextpow2(2*n-2); 118 else 119 nfft = 2*n-2; 120 end 121 nf = nfft/2;% number of frequencies 122 123 ACF = [ACF;zeros(nfft-2*n+2,1);ACF(n-1:-1:2)]; 124 125 Rper = real(fft(ACF,nfft));% periodogram 126 127 k = find(Rper<0); 128 if any(k) 129 % disp('Warning: negative spectral estimates') 130 % disp('Apply the parzen windowfunction ') 131 % disp('to the ACF in order to avoid this.') 132 % disp('The returned result is now only an approximation.') 133 %min(Rper(k)) 134 Rper(k)=0; % truncating negative values to ensure that 135 % that high frequency noise is not added to 136 % the Spectrum 137 end 138 139 140 k = find(Rper/max(Rper(:))<trunc); 141 if any(k) 142 Rper(k)=0; % truncating small values to ensure that 143 % that high frequency noise is not added to 144 % the Spectrum 145 end 146 147 %size(nf),size(dT) 148 S = createspec(spectype,ftype); 149 S.tr = R.tr; 150 S.h = R.h; 151 S.norm = R.norm; 152 S.note = R.note; 153 %fn = linspace(0,1, 154 switch ftype 155 case {'w'} 156 S.w = [0:(nf)]'/nf*pi/dT; % (rad/s) 157 S.S = abs(Rper(1:(nf+1)))*dT/pi; % (m^2*s/rad) one-sided spectrum 158 case {'f'} % spectype == sf 159 S.f =[0:(nf)]'/nf/2/dT; % frequency Hz if dT is in seconds 160 S.S =2*abs(Rper(1:(nf+1)))*dT; % (m^2*s) one sided spectrum 161 case {'k'} 162 S.k =[0:(nf)]'/nf*pi/dT;% (rad/m) 163 S.S = abs(Rper(1:(nf+1)))*dT/pi; % (m^3/rad) one-sided wavenumber spectrum 164 end 165 166 if rate>1 167 w = getfield(S,ftype); 168 wi = linspace(w(1),w(end),nf*rate).'; 169 S.S = interp(w,S.S,wi); 170 S = setfield(S,ftype,wi); 171 end 172 173 174 175 176 177
Comments or corrections to the WAFO group