I have to implement from scratch on MATLAB the periodogram of a very simple stochastic process:
$$X(t)=A_1\cos(2\pi f_1t+\phi_1)+A_2\cos(2\pi f_2t+\phi_2)$$
where $f_1=5$, $f_2=10$, $A_1,A_2\sim Rayleigh(\sigma^2)$, $\sigma^2=2$ and $\phi_1,\phi_2\sim U[0,2\pi]$.
To simulate it, I built a function that takes as inputs the frequencies $f_1,f_2,\sigma^2$, and allows for the user to define how many cycles of the periodic process they want to see, and how many data points. As this is a sum of two cosines, the fundamental frequency is $\min\{f_1,f_2\}$.
Therefore, my time vector is created as time=linspace(0,cycles/fundamental,datapoints)'.
I should then create a function that implements the periodogram of this stochastic process. It mainly consists of two steps:
Fourier transform of the data:
$$\mathcal{X}(f)=\sum_{t=0}^{n-1}x(t)e^{-2\pi ft}$$
Periodogram, as follows:
$$\hat{R}_x(f)=\frac{1}{n}|\mathcal{X}(f)|^2$$
I can say that my stochastic process was sampled with period Ts=(cycles/fundamental)/datapoints.
If I did everything right so far, my big doubt is: how should I define the frequency range?
On MATLAB's reference for FFT, the frequency range is: Fs*(0:datapoints/2)/datapoints, which seems reasonable, but in my case, if I do that, I may cause my spectrum "spikes" to occur not exactly at $f_1$ and $f_2$ unless my datapoints get very high.
On the other hand, pushing my datapoints to a large number causes my spectrum to show a huge range of frequencies that carry no relevant information at all, since my process is composed only by two low frequency cosines.
I believe that once I get this frequency range issue solved, the rest is straightforward:
w=-1i*2*pi*freq;
wt=w.*t;
W=exp(wt);
X_f=x'*W;
R_f=(1/(length(x)).*abs(X_f).^2;
Hope you guys can enlighten me with this issue.
By the way, the reference I am using is Georg Lindgren's Stationary stochastic Processes for Scientists and Engineers