# Spectral Analysis with fft

5 views (last 30 days)
Garfield17 on 17 Feb 2021
Answered: Jaynik on 19 Jul 2024
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

Jaynik on 19 Jul 2024
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);
figure
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
figure
hold on
legend('theoretical','fft')
Hope this help!

### Categories

Find more on Signal Generation and Preprocessing in Help Center and File Exchange

R2018b

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!