spectra.m (1120B)
1 % Spectra.m 2 % Computes an estimate of the spectral density using the 3 % Daniell filter to smooth the periodogram. 4 % 5 % h is normalized around one. 6 % pi/m = width of the window, e.g. m=16 7 8 function h = spectra(x,M); 9 10 x=demean(x); 11 y = fft(x); 12 Ix = y.*conj(y); 13 TF = length(Ix); 14 Ix = Ix/TF; 15 16 % Now apply the Daniell Filter 17 h=zeros(TF,1); 18 19 for t=1:TF; 20 hh=0; 21 dlta=floor(TF/(2*M)); 22 for i=1:(2*dlta+1); 23 j=i-1-dlta; % j is offset for smoothing 24 if t+j<=0; j=j+TF; end; % Since period = 2*pi 25 if t+j>TF; j=j-TF; end; 26 hh=hh+Ix(t+j); 27 end; % i 28 h(t) = hh/(2*dlta+1); 29 end % t 30 31 % Let's return the normalized spectral density: i.e. divide by variance 32 h=h(1:TF/2); 33 h=h/((x'*x)/length(x)); 34 35 % Plot a graph 36 freq=(1:TF/2)/TF; 37 plot(freq,h); 38 title('Spectral Density'); 39 40 % Finally, we use some specialized features of MATLAB's 41 % "Handle Graphics" to change the labeling on the x-axis 42 % in the bottom graph from frequency to its reciprocal. 43 44 % f = get(gca,'xtick'); 45 f=[0 .1 .2 .3 .4 .5]; 46 c = zeros(length(f),5); 47 for i = 1:length(f) 48 c(i,1:5) = sprintf('%5.1f',1/f(i)); 49 end 50 set(gca,'xticklabels',c) 51 xlabel('years/cycle')