Spectral Analysis with fft

Question

Garfield17 2021-2-17

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/747702-spectral-analysis-with-fft

回答： Jaynik 2024-7-19，13:45

Hi guys,

currently, I am trying to understand how spectral analysis of time series can be done in matlab, but unfortunately I have some difficulties in understanding the concept.

As a first stept I would like to compare the theoretical spectral density of an AR(1)-process with the one implied by a realization of an AR(1)-process.

This is what I have so far

rng(786786)
AR1 = arima('Constant',0,'AR',{0.9},'Variance',0.01^2);
Y1 = simulate(AR1,1000);
ts=Y1(end-200+1:end); %Take last 200 observations
T = length(ts); %Length of time series
%Calculate theoretical spectral function 
k=1:1:T-1;
theo_gamma=[0.01^2/(1-0.9^2) (0.9.^k)*0.01^2/(1-0.9^2)]'; %theoretical autocovariance fuction of an AR(1)-process
omega=2*pi*k'./T;
spec_theo=zeros(1,length(k));
for j=k
    spec_theo(j)=theo_gamma(1)/(2*pi)+1/pi*theo_gamma(2:end)'*cos(omega(j)*k');
end
%Calculating spectral density with fft
fft_ts = fft(ts);
spec_fft = abs(fft_ts).^2/(2*pi*T);
N=floor(T/2)+1;
x=1:N;
figure
plot(x,spec_theo(1:N))
hold on
plot(x,spec_fft(2:N+1)) %cutoff the first value
legend('theoretical','fft')

I divided the the period length of 2*pi into equidistant subintervals according to the time series length. My first question: Is the code correct?

Further, I would like to analyze some real data, so imagine that the 200 Realizations of the AR(1)-process correspond to 50 Years of quaterly data. In some textbooks, the are graphs showing the units of the x-Axis as "Cycles per year" or "Periods per Cycle". What do I have to do in order to get such graphs in my example, if we interpret the data as quaterly?

Thanks a lot!

Frank

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Jaynik 2024-7-19，13:45

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/747702-spectral-analysis-with-fft#answer_1488071

在 MATLAB Online 中打开

Hi Frank,

The code seems correct as you have correctly implemented the theoretical autocovariance function of an AR(1) process and used it to calculate the theoretical spectral function. You have also correctly implemented the Fast Fourier Transform to calculate the spectral dsitance of time series.

If you want to interpret the data as quaterly and want the units of the x-axis as "Cycles per year" or "Periods per Cycle", you would need to adjust the frequency accordingly.

For "cycles per year", you can divide the frequencies by the 4. Then we can adjust the x-axis values to match the quaterly data. Here is the code for plotting this:

freq_cycles_per_year = (0:N-1) / (T/4);
x_adj = x/4; % Adjust x for quarterly data
figure
plot(x_adj, spec_theo(1:N))
hold on
plot(x_adj, spec_fft(2:N+1)) %cutoff the first value
legend('theoretical','fft')

Similarly, for "periods per cycle", we can take the reciprocal of the frequencies to get periods. Here is the code for the same:

periods_per_cycle = 1 ./ freq_cycles_per_year;
periods_per_cycle(1) = NaN; % Avoid division by zero for the first element
x_adj = 4./x; % Adjust x for quarterly data
figure
plot(x_adj,spec_theo(1:N))
hold on
plot(x_adj,spec_fft(2:N+1)) %cutoff the first value
legend('theoretical','fft')