periodog.m

来自「Modeling and Forecasting Electricity Loa」· M 代码 · 共 46 行

function [power,freq,maxp]=periodog(x,fontsize)
%PERIODOG Periodogram of a time series.
%   [POWER,FREQ,MAXP]=PERIODOG(X) returns the POWER and FREQuency such that 
%   PLOT(FREQ,POWER) is the periodogram (estimate of the spectral density). 
%   Period (i.e. 1/FREQ) MAXP of the maximum value of the periodogram is 
%   also returned.
%   If no output parameters are supplied the periodogram is plotted 
%   automatically. PERIODOG(X,FONTSIZE) allows to specify a fontsize 
%   different than 14.

%   Written by Rafal Weron (2000.12.27)
%   Copyright (c) 2000-2006 by Rafal Weron

if nargin<2, fontsize = 14; end;

if length(x)/2 ~= floor(length(x)/2), 
   x = x(1:end-1);
end;

Y = fft(x);

N = length(Y);
Y(1) = [];
power = (abs(Y(1:N/2)).^2)/length(x);
nyquist = 1/2;
freq = ((1:N/2)/(N/2)*nyquist)';
period = 1./freq;
[max_p ind_p] = max(power);
maxp = period(ind_p);

if nargout<1,
   subplot(2,1,1)
   plot(freq,power,'.-')
   grid on
   set(gca,'Box','on','fontsize',fontsize);
   xlabel('Frequency','fontsize',fontsize);
   ylabel('Power','fontsize',fontsize);
   subplot(2,1,2)
   plot(1./freq,power,'.-')
   grid on
   set(gca,'Box','on','fontsize',fontsize);
   xlabel('Period','fontsize',fontsize);
   ylabel('Power','fontsize',fontsize);
end;