%WAVE_SIGNIF  Significance testing for the 1D Wavelet transform WAVELET%%   [SIGNIF,FFT_THEOR] = ...%      wave_signif(Y,DT,SCALE,SIGTEST,LAG1,SIGLVL,DOF,MOTHER,PARAM)%% INPUTS:%%    Y = the time series, or, the VARIANCE of the time series.%        (If this is a single number, it is assumed to be the variance...)%    DT = amount of time between each Y value, i.e. the sampling time.%    SCALE = the vector of scale indices, from previous call to WAVELET.%%% OUTPUTS:%%    SIGNIF = significance levels as a function of SCALE%    FFT_THEOR = output theoretical red-noise spectrum as fn of PERIOD%%% OPTIONAL INPUTS:% *** Note *** setting any of the following to -1 will cause the default%               value to be used.%%    SIGTEST = 0, 1, or 2.    If omitted, then assume 0.%%         If 0 (the default), then just do a regular chi-square test,%             i.e. Eqn (18) from Torrence & Compo.%         If 1, then do a "time-average" test, i.e. Eqn (23).%             In this case, DOF should be set to NA, the number%             of local wavelet spectra that were averaged together.%             For the Global Wavelet Spectrum, this would be NA=N,%             where N is the number of points in your time series.%         If 2, then do a "scale-average" test, i.e. Eqns (25)-(28).%             In this case, DOF should be set to a%             two-element vector [S1,S2], which gives the scale%             range that was averaged together.%             e.g. if one scale-averaged scales between 2 and 8,%             then DOF=[2,8].%%    LAG1 = LAG 1 Autocorrelation, used for SIGNIF levels. Default is 0.0%%    SIGLVL = significance level to use. Default is 0.95%%    DOF = degrees-of-freedom for signif test.%         IF SIGTEST=0, then (automatically) DOF = 2 (or 1 for MOTHER='DOG')%         IF SIGTEST=1, then DOF = NA, the number of times averaged together.%         IF SIGTEST=2, then DOF = [S1,S2], the range of scales averaged.%%       Note: IF SIGTEST=1, then DOF can be a vector (same length as SCALEs),%            in which case NA is assumed to vary with SCALE.%            This allows one to average different numbers of times%            together at different scales, or to take into account%            things like the Cone of Influence.%            See discussion following Eqn (23) in Torrence & Compo.%%%----------------------------------------------------------------------------%   Copyright (C) 1995-1998, Christopher Torrence and Gilbert P. Compo%   University of Colorado, Program in Atmospheric and Oceanic Sciences.%   This software may be used, copied, or redistributed as long as it is not%   sold and this copyright notice is reproduced on each copy made.  This%   routine is provided as is without any express or implied warranties%   whatsoever.%----------------------------------------------------------------------------function [signif,fft_theor] = ...	wave_signif(Y,dt,scale1,sigtest,lag1,siglvl,dof,mother,param);if (nargin < 9), param = -1;, endif (nargin < 8), mother = -1;, endif (nargin < 7), dof = -1;, endif (nargin < 6), siglvl = -1;, endif (nargin < 5), lag1 = -1;, endif (nargin < 4), sigtest = -1;, endif (nargin < 3)	error('Must input a vector Y, sampling time DT, and SCALE vector')endn1 = length(Y);J1 = length(scale1) - 1;scale(1:J1+1) = scale1;s0 = min(scale);dj = log(scale(2)/scale(1))/log(2.);if (n1 == 1)	variance = Y;else	variance = std(Y)^2;endif (sigtest == -1), sigtest = 0;, endif (lag1 == -1), lag1 = 0.0;, endif (siglvl == -1), siglvl = 0.95;, endif (mother == -1), mother = 'MORLET';, endmother = upper(mother);% get the appropriate parameters [see Table(2)]if (strcmp(mother,'MORLET'))  %----------------------------------  Morlet	if (param == -1), param = 6.;, end	k0 = param;	fourier_factor = (4*pi)/(k0 + sqrt(2 + k0^2)); % Scale-->Fourier [Sec.3h]	empir = [2.,-1,-1,-1];	if (k0 == 6), empir(2:4)=[0.776,2.32,0.60];, endelseif (strcmp(mother,'PAUL'))  %--------------------------------  Paul	if (param == -1), param = 4.;, end	m = param;	fourier_factor = 4*pi/(2*m+1);	empir = [2.,-1,-1,-1];	if (m == 4), empir(2:4)=[1.132,1.17,1.5];, endelseif (strcmp(mother,'DOG'))  %---------------------------------  DOG	if (param == -1), param = 2.;, end	m = param;	fourier_factor = 2*pi*sqrt(2./(2*m+1));	empir = [1.,-1,-1,-1];	if (m == 2), empir(2:4) = [3.541,1.43,1.4];, end	if (m == 6), empir(2:4) = [1.966,1.37,0.97];, endelse	error('Mother must be one of MORLET,PAUL,DOG')endperiod = scale.*fourier_factor;dofmin = empir(1);     % Degrees of freedom with no smoothingCdelta = empir(2);     % reconstruction factorgamma_fac = empir(3);  % time-decorrelation factordj0 = empir(4);        % scale-decorrelation factorfreq = dt ./ period;   % normalized frequencyfft_theor = (1-lag1^2) ./ (1-2*lag1*cos(freq*2*pi)+lag1^2);  % [Eqn(16)]fft_theor = variance*fft_theor;  % include time-series variancesignif = fft_theor;if (dof == -1), dof = dofmin;, endif (sigtest == 0)    % no smoothing, DOF=dofmin [Sec.4]	dof = dofmin;	chisquare = chisquare_inv(siglvl,dof)/dof;	signif = fft_theor*chisquare ;  % [Eqn(18)]elseif (sigtest == 1)  % time-averaged significance	if (length(dof) == 1), dof=zeros(1,J1+1)+dof;, end	truncate = find(dof < 1);	dof(truncate) = ones(size(truncate));	dof = dofmin*sqrt(1 + (dof*dt/gamma_fac ./ scale).^2 );   % [Eqn(23)]	truncate = find(dof < dofmin);	dof(truncate) = dofmin*ones(size(truncate));   % minimum DOF is dofmin	for a1 = 1:J1+1		chisquare = chisquare_inv(siglvl,dof(a1))/dof(a1);		signif(a1) = fft_theor(a1)*chisquare;	endelseif (sigtest == 2)  % time-averaged significance	if (length(dof) ~= 2)		error('DOF must be set to [S1,S2], the range of scale-averages')	end	if (Cdelta == -1)		error(['Cdelta & dj0 not defined for ',mother, ...			' with param = ',num2str(param)])	end	s1 = dof(1);	s2 = dof(2);	avg = find((scale >= s1) & (scale <= s2));  % scales between S1 & S2	navg = length(avg);	if (navg == 0)		error(['No valid scales between ',num2str(s1),' and ',num2str(s2)])	end	Savg = 1./sum(1 ./ scale(avg));       % [Eqn(25)]	Smid = exp((log(s1)+log(s2))/2.);     % power-of-two midpoint	dof = (dofmin*navg*Savg/Smid)*sqrt(1 + (navg*dj/dj0)^2);  % [Eqn(28)]	fft_theor = Savg*sum(fft_theor(avg) ./ scale(avg));  % [Eqn(27)]	chisquare = chisquare_inv(siglvl,dof)/dof;	signif = (dj*dt/Cdelta/Savg)*fft_theor*chisquare;    % [Eqn(26)]else	error('sigtest must be either 0, 1, or 2')endreturn% end of code