Home > wafo > onedim > dat2spec2.m

dat2spec2

PURPOSE ^

Estimate one-sided spectral density, version 2.

SYNOPSIS ^

S = dat2spec2(xn,L,g,plotflag,wdef,dT,chopOffHighFreq)

DESCRIPTION ^

 DAT2SPEC2 Estimate one-sided spectral density, version 2.  
       using a Parzen window function on  
       the estimated autocovariance function. 
   
  CALL:  S = dat2spec2(x,L,g,plotflag,wdef,dT) 
  
          S = A structure containing: 
              S    = spectral density 
              w    = angular frequency 
              tr   = transformation g 
              h    = water depth (default inf) 
              type = 'freq' 
              note = Memorandum string 
              date = Date and time of creation  
              L    = maximum lag size of the window function.  
              Bw   = Bandwidth of the smoothing window which is used  
                     in the estimated spectrum. (rad/sec or Hz) 
              
  
          x =  m column data matrix with sampled times in the first column 
               and values the next columns.     
  
          L = maximum lag size of the window function.  
              If no value is given, the lag size is set to 
              be the lag where the auto correlation is less than  
              2 standard deviations. (maximum 300)  
                           
          g = the transformation assuming that x is a sample of a  
              transformed Gaussian process. If g is empty then 
              x  is a sample of a Gaussian process (Default) 
  
   plotflag = 1 plots the spectrum, S,  
              2 plot 10log10(S) and 
              3 plots both the above plots 
  
         wdef = 'parzen'   then a Parzen window is used (default). 
                'hanning'  then a Hanning window is used.  
  
          dT  = sampling interval. Default is dT=x(2,1)-x(1,1) or 1Hz. 
  
   As L decreases the estimate becomes smoother and Bw increases.  
   If we want to resolve peaks in S which is Bf (Hz or rad/sec) apart  
   then Bw < Bf.  
    
  See also   dat2tr, dat2cov, dat2spec

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

001 function S = dat2spec2(xn,L,g,plotflag,wdef,dT,chopOffHighFreq) 
002 %DAT2SPEC2 Estimate one-sided spectral density, version 2.  
003 %      using a Parzen window function on  
004 %      the estimated autocovariance function. 
005 %  
006 % CALL:  S = dat2spec2(x,L,g,plotflag,wdef,dT) 
007 % 
008 %         S = A structure containing: 
009 %             S    = spectral density 
010 %             w    = angular frequency 
011 %             tr   = transformation g 
012 %             h    = water depth (default inf) 
013 %             type = 'freq' 
014 %             note = Memorandum string 
015 %             date = Date and time of creation  
016 %             L    = maximum lag size of the window function.  
017 %             Bw   = Bandwidth of the smoothing window which is used  
018 %                    in the estimated spectrum. (rad/sec or Hz) 
019 %             
020 % 
021 %         x =  m column data matrix with sampled times in the first column 
022 %              and values the next columns.     
023 % 
024 %         L = maximum lag size of the window function.  
025 %             If no value is given, the lag size is set to 
026 %             be the lag where the auto correlation is less than  
027 %             2 standard deviations. (maximum 300)  
028 %                          
029 %         g = the transformation assuming that x is a sample of a  
030 %             transformed Gaussian process. If g is empty then 
031 %             x  is a sample of a Gaussian process (Default) 
032 % 
033 %  plotflag = 1 plots the spectrum, S,  
034 %             2 plot 10log10(S) and 
035 %             3 plots both the above plots 
036 % 
037 %        wdef = 'parzen'   then a Parzen window is used (default). 
038 %               'hanning'  then a Hanning window is used.  
039 % 
040 %         dT  = sampling interval. Default is dT=x(2,1)-x(1,1) or 1Hz. 
041 % 
042 %  As L decreases the estimate becomes smoother and Bw increases.  
043 %  If we want to resolve peaks in S which is Bf (Hz or rad/sec) apart  
044 %  then Bw < Bf.  
045 %   
046 % See also   dat2tr, dat2cov, dat2spec 
047  
048 % May chop off high frequencies in order to  
049 % get the same irregularity factor in the spectrum as  
050 % in x. 
051  
052  
053 % References: 
054 % Georg Lindgren and Holger Rootzen (1986) 
055 % "Stationära stokastiska processer",  pp 173--176. 
056 %  
057 % Gareth Janacek and Louise Swift (1993) 
058 % "TIME SERIES forecasting, simulation, applications", 
059 % pp 75--76 and 261--268 
060 % 
061 % Emanuel Parzen (1962), 
062 % "Stochastic Processes", HOLDEN-DAY, 
063 % pp 66--103 
064  
065  
066 % Tested on: Matlab 5.3 
067 % history: 
068 % Modified by Per A. Brodtkorb 14.08.98,25.05.98 
069 %  - add a nugget effect to ensure that round off errors 
070 %    do not result in negative spectral estimates 
071 % Modified by svi 29.09.99 
072 %  - The program is not estimating transformation g any more. 
073 % modified pab 22.10.1999 
074 %  - fixed so that x in fact can be a m column matrix 
075 %  - updated info put into the spectral structure. 
076 %  - updated help header 
077 % modified pab 03.11.1999 
078 %  - fixed a bug: line 152 wrong array dim. when m>2 
079 % modified jr 29.01.2000 
080 %  - modified from dat2spec, new name: dat2spec2 
081 %  - lines related to confidence bands removed 
082 %  - call to psd removed 
083 % Revised pab Dec 2003 
084 % -fixed a bug for wdef and updated helpheader accordingly. 
085 %   
086  
087 nugget=10^-12; 
088  
089 freqtype='w'; %options are 'f' and 'w' 
090  
091 xx=xn; 
092 istime=1; 
093 [n m]= size(xx); 
094  
095 if min(m,n)==1,  
096   xx=[ (1:n)' xx(:)];istime=0; 
097   n=max(m,n); 
098   m=2; 
099 end 
100  
101 if (nargin < 7) | isempty(chopOffHighFreq), 
102   chopOffHighFreq=0;  
103   % chop off high frequencies in order to get the same  
104   % irregularity factor in the spectrum as in the data 
105   % not a good idea 
106 end 
107  
108 if (nargin < 6) | isempty(dT), 
109   dT=xx(2,1)-xx(1,1); 
110   if ~istime 
111     disp(' Warning: The sampling frequency is undetermined and set to 1 Hz.') 
112   end 
113 end 
114  
115 if (nargin < 5) | isempty(wdef) 
116   wdef=1; % 1=parzen window 2=hanning window 
117 elseif ischar(wdef) 
118   switch lower(wdef(1)) 
119    case 'p', 
120     wdef = 1; % parzen 
121    case 'h', 
122     wdef = 2; % hanning 
123    otherwise 
124     error('unknown option') 
125   end 
126 end 
127  
128 if (nargin < 4) | isempty(plotflag) 
129   plotflag = 0;  
130 end 
131 if (nargout == 0) & (plotflag==0) 
132   plotflag = 1;  
133 end 
134  
135 if (nargin<3)  
136   g=[]; 
137 end 
138  
139 if  (isempty(g)),  
140   yy=xx; 
141 else 
142   yy = dat2gaus(xx,g); 
143 end 
144  
145 % By using a tapered data window to smooth the data at each 
146 % end of the record has the effect of sharpening 
147 % the spectral window. 
148 % NB! the resulting spectral estimate must be  
149 % normalized in order to correct for the loss of 
150 % amplitude (energy) caused by the data taper. 
151 %taper=bingham(n); 
152 %yy(:,2)=bingham(n).*yy(:,2); 
153 %display('****1'); 
154 %break; 
155 yy(:,2:m)=(yy(:,2:m)-ones(n,1)*mean(yy(:,2:m)) );%/std(yy(:,2)); 
156  
157 max_L=300; 
158  
159 changed_L =0; 
160 if (nargin < 2 | isempty(L)) 
161   L=min(n-2,ceil(4/3*max_L)); 
162 %else 
163   changed_L =1; 
164 end 
165  
166 %if ~strcmp(method,'psd') | changed_L,  
167   r=0; 
168   stdev=0; 
169   for ix=2:m 
170     R=dat2cov(yy(:,[1 ix])); 
171     r=r+R.R(:); 
172     stdev=stdev+R.stdev(:); 
173   end 
174   R.stdev=stdev/(m-1); 
175   r=r/(m-1); 
176   R.R=r; 
177   if   changed_L, 
178     % finding where ACF is less than 2 st. deviations. 
179     L=find(abs(r(1:max_L))>2*R.stdev(1:max_L))+1; % a better L value   
180     L=min(L(end),max_L); 
181   end   
182 %end 
183  
184 if wdef==1             % Parzen window 
185   if changed_L,        % modify L so that hanning and Parzen give appr. the same  
186                        % result 
187     L=min(floor(4*L/3),n-2);disp(['The default L is set to ' num2str(L) ]) 
188   end 
189   v=floor(3.71*n/L);   % degrees of freedom used in chi^2 distribution 
190   win=parzen(2*L-1);   % Does not give negative estimates of the spectral density 
191   Be=2*pi*1.33/(L*dT); % bandwidth (rad/sec) 
192    
193 elseif wdef==2         % Hanning window 
194   if changed_L, 
195     disp([' The default L is set to ' num2str(L) ]) 
196   end 
197   v=floor(2.67*n/L);   % degrees of freedom used in chi^2 distribution 
198   win=hanning(2*L-1);  % May give negative estimates of the spectral density 
199   Be =2*pi/(L*dT);     % bandwidth (rad/sec) 
200 end 
201  
202 nf=2^nextpow2(2*L-2);  %  Interpolate the spectrum with rate = 2 
203 nfft=2*nf; 
204  
205 S=createspec('freq',freqtype); 
206 S.tr=g; 
207 S.note=['dat2spec(',inputname(1),')']; 
208 S.L=L; 
209 S.S=zeros(nf+1,m-1); 
210  
211 % add a nugget effect to ensure that round off errors 
212 % do not result in negative spectral estimates 
213  
214 r=r+nugget; 
215 rwin=r(1:L).*win(L:(2*L-1));  
216 Rper=real(fft([rwin; zeros(nfft-(2*L-1),1); rwin(L:-1:2)],nfft)); 
217  
218 rpmin=min(Rper); 
219 if rpmin<0 
220   disp(' Warning: negative spectral estimates') 
221 end 
222  
223 if strcmp(freqtype,'w') 
224   S.w=[0:(nf)]'/nf*pi/dT;              % (rad/s) 
225   S.S=abs(Rper(1:(nf+1),1))*dT/pi;     % (m^2*s/rad) one-sided spectrum 
226   S.Bw=Be; 
227 else % freqtype == f 
228   S.f=[0:(nf)]'/nf/2/dT;               % frequency Hz if dT is in seconds 
229   S.S=2*abs(Rper(1:(nf+1),1))*dT;      % (m^2*s) one-sided spectrum 
230   S.Bw=Be/(2*pi);                      % bandwidth in Hz 
231 end 
232  
233  
234   
235  
236 N=floor(nf/10); 
237 % cutting off high frequencies  
238 % in this way may not be a very good idea. 
239  
240 if (N>3)&chopOffHighFreq,  
241    
242     % The data must be Gaussian in order for this proc to be correct. 
243     [g0 test cmax irr]  = dat2tr(xx,'nonlinear'); 
244   ind=(nf-4):(nf+1); 
245   [S,m4,m2]=wnormspec(S,0); 
246  
247   while (sqrt(m4/m2^2) > irr) & (ind(1)>1)  
248     S.S(ind)=0; 
249     [S,m4,m2]=wnormspec(S,0); 
250     ind=ind-5; 
251   end 
252   fcut=S.w(min(ind(end)+5,nf)); % cut off frequency 
253   if ind(1)<1, 
254     disp('DAT2SPEC: Error in cutting off high frequencies, try other L-values') 
255   end 
256 end 
257  
258  
259  
260 %-----------------------------------  
261 %  
262 %     Plotting the Spectral Density  
263 % 
264 %----------------------------------- 
265  
266 if plotflag>0 
267  wspecplot(S,plotflag) 
268 end 
269  
270

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