Home > wafo > spec > cov2spec.m

cov2spec

PURPOSE ^

Computes spectral density given the auto covariance function

SYNOPSIS ^

S = cov2spec(R,rate,nugget,trunc,fast);

DESCRIPTION ^

 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

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

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

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