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
Normal probability density function | |
Control axis scaling and appearance. | |
Display message and abort function. | |
Get handle to current axis. | |
Get object properties. | |
Linear plot. | |
Fit polynomial to data. | |
Set object properties. | |
Stairstep plot. | |
Graph title. | |
X-axis label. | |
Y-axis label. |
Calculates the number of upcrossings from a cycle count | |
Extracts level-crossing spectrum from min2Max cycles. | |
Quick test of the routines in module 'cycles' |
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
Comments or corrections to the WAFO group