Stack Exchange
FFT of real valued even vector in MATLAB is not real
I am trying to find the PSD of a signal by taking the DFT of its auto-covariance sequence.
My code is as follows
% define constants
m = 256;
k = 0:m - 1;
a1 = -0.9;
b1 = 0.8;
sigma = 1;

% compute ACS for positive lags k
r = (sigma^2 / (1 - abs(a1)^2)) * ( (1 + abs(b1)^2) * d(k) + b1 * d(k - 1) + conj(b1) * d(k + 1));
% compute acs from -(m-1) to +(m-1)
acs = [flip(r(2:end)) r];

% compute PSD
phi = fft(acs);


% auxilliary function that helps in the computation of the ACS
% it computes the function
% d(k) = (-a1)^k, k >=0,
%        (-a1*)^(-k), k < 0
function d = d(x)
    a1 = -0.9;
    x1 = x >= 0;
    x2 = ~x1;

    d = x1 .* (-a1).^x + x2 .* (-conj(a1)).^(-x);
end

The acs variable is plotted as

and can be seen to be even. However when i take the fft of the acs the results i get are:
for the real part

and for the imaginary part


Probably a discretisation or numerical issue. Your real part is on the order of 300, whereas your numeric part is only on the order of 8. That's a pretty large difference.

Shouldn't the real part also be positive since by definition the PSD is always positive?

You mean that I have to move acs(256), which is basically r[0], to acs(1),  followed by r[k] for positive lags and then r[k] for negative lags?

Only the DC and alternating components of a real to complex conjugate symmetric Fourier transform are guaranteed to be real! If you want the FFT of your presumed symmetric function to be real you either have to centre it on the origin of the array (for standard FFT) or use an FFT that puts the centre component at N/2. All classical  real to complex conjugate symmetric FFTs assume  their origin is at x[0]. The power spectrum you requires you to compute the magnitude of each complex component. While you have phase info you can reconstruct the original waveform exactly.

Try acs = [r, flip(r(2:end))];

Cross-posted at https://dsp.stackexchange.com/q/89821/33605 — please don’t cross-post.

Already has an answer here.
