%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % Spectral analysis of a contrived time series % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Start by constructing a time series we'd like to evaluate % We are designing these so x has unit variance. ntim = 1000; % The time series will have 1500 points % Use one of the following two routines - generate an AR1 OR an AR2 % process. % AR1: if false; x = zeros(ntim,1); % Initialize x x(1) = randn(1); % Initial condition - random rho1 = 0.8; % One lag autocorrelation for i = 2:ntim; x(i) = rho1*x(i-1)+sqrt(1-rho1^2)*randn(1); end end % AR2: if true; x = zeros(ntim,1); % Initialize x x(1) = randn(1); % Initial condition - random rho1 = 0.5; % One lag autocorrelation rho2 = -0.35; % Two lag autocorrelation % Get parameters alpha1 and alpha2 via Yule-Walker equations fi = [rho1; rho2]; Big_Fi = [1 rho1; rho1 1]; % X has unit variance alph = Big_Fi \ fi; % Solve Yule-Walker equations alpha1 = alph(1); alpha2 = alph(2); % Now, generate time series noiseamp = sqrt(1 - alpha1*rho1 - alpha2*rho2); for i = 3:ntim; x(i) = alpha1*x(i-1)+alpha2*x(i-2)+noiseamp*randn(1); end end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Let's analyze this using overlapping windows %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% nfft = 128; % Length of window noverlap = nfft/2; % Amount of overlap per segment %window = boxcar(nfft); % Window we'll use ... window = hanning(nfft); % Or this one %window = hamming(nfft); % Or this one Fs = 1/1; % Sampling rate: observations per day (in the case of pentad % data, we might use the following): % Fs = 0.2 = (1 observation)/(5 days); % Or, for monthly data: % Fs = 12 = (12 observations)/(1 year); DOF = 2*ntim/nfft; % Degrees of Freedom per spectral estimate % Calculate the spectrum [p,f] = spectrum(x, nfft, noverlap, window, Fs); % Generate null hypotheses: h = [0:(nfft/2)]; % AR1: if false; pAR = (1-rho1^2)./(1+rho1^2-2*rho1*cos(2*pi*h/nfft)); end % AR2: if true; g = alpha1*(1-alpha2)*cos(2*pi*h/nfft)+alpha2*cos(4*pi*h/nfft); pAR = noiseamp^2./(1+alpha1^2+alpha2^2-2*g); end % Normalize to make sure everything has equal variance (noting that the % first two spectral estimates are typically strongly affected by the % window). pAR = pAR * sum(p(3:end,1))/sum(pAR(3:end)); % Get f-value: fsig = finv(0.95, DOF, 100000); % If you do not have the stats toolbox, % then just look this up in a table. % Plot the spectrum in a number of ways: figure(1); orient landscape; clf; % Linear subplot(2,2,1); h = plot(f, p(:,1), '-k', f, pAR, '-r', f, pAR*fsig, '--r'); set(h(1), 'linewidth', 2); % Make the spectrum a thicker line axis([0 0.5 -.1 14.1]); set(gca, 'XTick', [0:.1:0.5], 'YTick', 0:2:14); grid on; t1 = title('Spectrum: Linear plot', 'fontsize', 12); xl = xlabel('Frequency (cycles per day)', 'fontsize', 10); yl = ylabel('Power', 'fontsize', 10); % Log-Linear subplot(2,2,2); h = semilogy(f, p(:,1), '-k', f, pAR, '-r', f, pAR*fsig, '--r'); set(h(1), 'linewidth', 2); % Make the spectrum a thicker line axis([0 0.5 -.1 14.1]); set(gca, 'XTick', [0:.1:0.5]); grid on; t1 = title('Spectrum: Log-linear plot', 'fontsize', 12); xl = xlabel('Frequency (cycles per day)', 'fontsize', 10); yl = ylabel('log_1_0(Power)', 'fontsize', 10); % Log-log subplot(2,2,3); h = loglog(f, p(:,1), '-k', f, pAR, '-r', f, pAR*fsig, '--r'); set(h(1), 'linewidth', 2); % Make the spectrum a thicker line axis([0 0.5 -.1 14.1]); grid on; set(gca, 'XTick', 2.^[-6:-1], 'XTickLabel', 2.^[6:-1:1]); set(gca, 'XMinorGrid', 'off'); t1 = title('Spectrum: Log-log plot', 'fontsize', 12); xl = xlabel('Period (days per cycle)', 'fontsize', 10); yl = ylabel('log_1_0(Power)', 'fontsize', 10); % f*log(P) vs. log(f) subplot(2,2,4); h = plot(log(f), f.*log(p(:,1)), '-k', ... log(f), f.*log(pAR'), '-r', ... log(f), f.*log(pAR'*fsig), '--r'); set(h(1), 'linewidth', 2); % Make the spectrum a thicker line axis([log(f(2)) log(0.5) -1.75 0.5]); grid on; set(gca, 'XTick', log(2.^[-6:-1]), 'XTickLabel', 2.^[6:-1:1]); set(gca, 'XMinorGrid', 'off'); t1 = title('Spectrum: f*ln(P) vs ln(f)', 'fontsize', 12); xl = xlabel('Period f (days per cycle)', 'fontsize', 10); yl = ylabel('f*ln(P)', 'fontsize', 10);