Error in power spectral density computation

Question

I've written a program to compute the power spectral density of a signal. The problem I'm having is that the computed spectrum is quite a bit more noisy than it should be (comparing to the spectrum computed by Matlab), especially at low frequencies.

http://i.imgur.com/fGOGT.png

The red in this plot shows the spectrum computed by Matlab's cpsd command (0 overlap, NFFT=window=8192, 200Hz signal), and the green shows my own. There's a constant offset (in the log domain) between them, but that's okay for now. What I'm more concerned about is that my implementation is noisier, especially in the 0-3Hz range. The algorithm I've implemented is

1. Zero out a vector of length NFFT/2+1, for averaging.
2. Until there are no more samples:
   1. Read NFFT samples (0-fill the remainder if there are less than NFFT)
   2. Compute the 1-dimensional real NFFT
   3. Compute the magnitude of each complex result, and divide by NFFT
   4. Add this to the corresponding position in the averaging vector
3. Divide each sum by the twice the sampling frequency times the number of windows read

Sizing a window different than NFFT, and overlap, aren't included in the algorithm, but when compared to the Matlab output for window=NFFT and no overlap it shouldn't matter.

Can anyone speak as to why this results in low frequency instability?

Answer 1

Matlab splits your signal into 8 overlapping segments (by default). The difference I see between Matlab and your algorithm is that Matlab also multiplies each segment by a Hamming window. You should get pretty close results if you do that.