Home > wafo > onedim > lcplot.m

lcplot

PURPOSE ^

Plots level-crossing spectrum (lc)

SYNOPSIS ^

h = lcplot(lc,plotflag,ma,sa);

DESCRIPTION ^

 LCPLOT Plots level-crossing spectrum (lc) 
 
  CALL:  h = lcplot(lc,plotflag,ma,sa);
 
         h = handle of the graphics object 
        lc = two column matrix with levels and number of upcrossings,
             i.e., level-crossing spectrum (LCS).
  plotflag = 1  plots LCS
             2  plots LCS and the theoretical one for a Gaussian
                process (default) 
             11 plot level-crossing intensity
             12 plot level-crossing intensity and the theoretical one for
                a Gaussian process  
   
     ma,sa = mean and standard deviation of the process 
             (default ma = 0, sa estimated from lc)
 
  Example: 
    x = load('sea.dat');
   lc = dat2lc(x,0.2); 
   lcplot(lc) 
 
  See also  dat2lc, mm2lc, wnormpdf, polyfit

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

001 function h = lcplot(lc,plotflag,ma,sa);
002 %LCPLOT Plots level-crossing spectrum (lc) 
003 %
004 % CALL:  h = lcplot(lc,plotflag,ma,sa);
005 %
006 %        h = handle of the graphics object 
007 %       lc = two column matrix with levels and number of upcrossings,
008 %            i.e., level-crossing spectrum (LCS).
009 % plotflag = 1  plots LCS
010 %            2  plots LCS and the theoretical one for a Gaussian
011 %               process (default) 
012 %            11 plot level-crossing intensity
013 %            12 plot level-crossing intensity and the theoretical one for
014 %               a Gaussian process  
015 %  
016 %    ma,sa = mean and standard deviation of the process 
017 %            (default ma = 0, sa estimated from lc)
018 %
019 % Example: 
020 %   x = load('sea.dat');
021 %  lc = dat2lc(x,0.2); 
022 %  lcplot(lc) 
023 %
024 % See also  dat2lc, mm2lc, wnormpdf, polyfit
025 
026 % Tested on: Matlab 6.0
027 % History:
028 % revised pab Feb2004  
029 %  - added plotflag ==11, 12
030 % revised by jr 06.07.00. 
031 % - Division by 2 from 05.07.00 removed (line 42). 
032 % - sa -> sa^2 (line 66)
033 % revised by jr 05.07.00. Line 39: division by 2
034 % revised by pab 11.08.99
035 % - removed the plot procedure from the old mm2cross function  
036 % - added the option of overplotting of the theoretical one 
037 %    for a Gaussian process
038 
039 error(nargchk(1,4,nargin))
040 if nargin<2|isempty(plotflag)
041  plotflag=2;
042 end
043 [cmax icmax]=max(lc(:,2));% maximum crossing
044 if (nargin <4|isempty(sa))& mod(plotflag,10)==2, 
045   % if stdev of x unknown then estimate it
046   cros   = lc;
047   indz   = (cros(:,2)==0);
048   lncros = cros; 
049   logcros(indz,2)  = inf; % this is done to avoid a warning message
050   logcros(~indz,2) = -log(cros(~indz,2));
051   logcmin      = logcros(icmax,2);
052   logcros(:,2) = sqrt(2*abs(logcros(:,2)-logcmin)); %sqrt(2*logcros)
053   logcros(1:icmax,2) = 2*logcros(icmax,2)-logcros(1:icmax,2); 
054   p = polyfit(lc(10:end-9,1),logcros(10:end-9,2),1); %least square fit
055   
056   sa = 1/p(1);% estimated standard deviation of x
057   
058   %plot(lc(:,1),logcros(:,2)),hold on
059   %plot(lc(:,1),polyval(p,lc(:,1)),'r'),hold off, pause
060 end
061 
062 if (nargin<3|isempty(ma))& mod(plotflag,10)==2, 
063   ma=0; % the mean is apriori zero.
064   %ma=p(2); %, lc(icmax,1) % new estimated mean
065   %ma=lc(icmax,1); % maximum cr intensity should be at the mean
066   % but this not always the case
067 end
068 
069 if plotflag>2
070   lc(:,2) = lc(:,2)/cmax;
071   cmax = 1;
072   ylabtxt = 'Relative crossing intensity';
073 else
074   ylabtxt = 'Number of upcrossings';
075 end
076 
077 %  clf
078 hh = stairs(lc(:,1),lc(:,2));
079 axis([-1.05*max(abs(lc(:,1))) 1.05*max(abs(lc(:,1))) 0 1.05*cmax])
080 title('Crossing spectrum')
081 xlabel('level u')
082 ylabel(ylabtxt)
083 
084 % Make a special graph box.
085 set(gca,'Box','off','TickLength',0.01*[1 1],'TickDir','out')
086 
087 
088 % bar(d(:,1),d(:,2))
089 if (mod(plotflag,10)==2)
090   y = wnormpdf(lc(:,1),ma,sa.^2);  
091   y = cmax*y/y(icmax);%Normalization needed to overplot the crossingspectrum.
092   np = get(gca,'NextPlot');    
093   set(gca,'NextPlot','add')    
094   hh1 = plot(lc(:,1),y,'r--','LineWidth',1);% Plots density line over histogram.
095   set(gca,'NextPlot',np) 
096 end
097 
098 if nargout == 1 & plotflag==2
099   h = [hh; hh1];
100 elseif nargout==1
101   h=hh;
102 end  
103  
104 return
105   
106

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