Frequency computation and fast fourier transform in Matlab

Question

I have a question related to Fast Fourier transform. I want to calculate the phase and make FFT to draw power spectral density. However when I calculate the frequency f, there are some errors. This is my program code:

n = 1:32768;

T = 0.2*10^-9; % Sampling period

Fs = 1/T; % Sampling frequency

Fn = Fs/2; % Nyquist frequency

omega = 2*pi*200*10^6; % Carrier frequency

L = 32768; % % Length of signal

t = (0:L-1)*T; % Time vector

x_signal(n) = cos(omega*T*n + 0.1*randn(size(n))); % Additive phase noise (random)

y_signal(n) = sin(omega*T*n + 0.1*randn(size(n))); % Additive phase noise (random)

theta(n) = atan(y_signal(n)/x_signal(n));

f = (theta(n)-theta(n-1))/(2*pi)

Y = fft(f,t);

PSD = Y.*conj(Y); % Power Spectral Density

%Fv = linspace(0, 1, fix(L/2)+1)*Fn; % Frequency Vector

SleuthEye SleuthEye · Accepted Answer · 2017-03-13T03:04:22

As posted, you would get the error

error: subscript indices must be either positive integers less than 2^31 or logicals

which refers to the operation theta(n-1) when n=1 which results in an index of 0 (which is out of bounds since Matlab uses 1-based indexing). To avoid that could use a subset of indices in n:

f = (theta(n(2:end))-theta(n(1:end-1)))/(2*pi);

That said, if you are doing this to try to obtain an instantaneous measure of the frequency, then you will have a few more issues to deal with. The most trivial one is that you should also divide by T. Not as obvious is the fact that as given, theta is a scalar due to the use of the / operator (see Matlab's mrdivide) rather than the ./ operator which performs element-wise division. So a better expression would be:

theta(n) = atan(y_signal(n)./x_signal(n));

Now, the next problem you might notice is that you are actually losing some phase information since the result of atan is [-pi/2,pi/2] instead of the full [-pi,pi] range. To avoid this you should instead be using atan2:

theta(n) = atan2(y_signal(n), x_signal(n));

Even with this, you are likely to notice that the estimated frequency regularly has spikes whenever the phase jumps between near -pi and near pi. This can be avoided by computing the phase difference modulo 2*pi:

f = mod(theta(n(2:end))-theta(n(1:end-1)),2*pi)/(2*pi*T);

A final thing to note: when calling the fft, you should not be passing in a time variable (the input is implicitly assumed to be sampled at regular time intervals). You may however specify the desired length of the FFT. So, you would thus compute Y as follow:

Y = fft(f, L);

And you could then plot the resulting PSD using:

Fv = linspace(0, 1, fix(L/2)+1)*Fn; % Frequency Vector
plot(Fv, abs(PSD(1:L/2+1)));

Frequency computation and fast fourier transform in Matlab

1 Answers