Thread Subject:

fft: small documentation bug?

Hello Matlab users,

I have used the code example in the documentation for the fft function

http://www.mathworks.de/help/techdoc/ref/fft.html

for my own program. The code claims to show an amplitude spectrum. The signal components shown in the example indicate this, and it can be easily tested with various numerical examples. There seems to be one exeption, though: The first element of the spectrum, 'the DC component'. The code

Y = fft(y,NFFT)/L;
f = Fs/2*linspace(0,1,NFFT/2+1);
plot(f,2*abs(Y(1:NFFT/2+1)))

appears wrong to me with respect to this element, because it is the only one that must not be multiplied by 2. I have tested this by adding a DC component to test signals, which appears with twice the amplitude in the spectrum. I would like to know what the experts say, because it could also be a matter of the definition of an amplitude spectrum, etc., etc.

Heinrich

"Heinrich Acker" wrote in message <jg950d$g5k$1@newscl01ah.mathworks.com>...
> Hello Matlab users,
>
> I have used the code example in the documentation for the fft function
>
> http://www.mathworks.de/help/techdoc/ref/fft.html
>
> for my own program. The code claims to show an amplitude spectrum. The signal components shown in the example indicate this, and it can be easily tested with various numerical examples. There seems to be one exeption, though: The first element of the spectrum, 'the DC component'. The code
>
> Y = fft(y,NFFT)/L;
> f = Fs/2*linspace(0,1,NFFT/2+1);
> plot(f,2*abs(Y(1:NFFT/2+1)))
>
> appears wrong to me with respect to this element, because it is the only one that must not be multiplied by 2. I have tested this by adding a DC component to test signals, which appears with twice the amplitude in the spectrum. I would like to know what the experts say, because it could also be a matter of the definition of an amplitude spectrum, etc., etc.
>
> Heinrich

You're right Heinrich, but even a bit more than that, the Nyquist only occurs once as well and therefore should not be doubled.

spectrum.periodogram handles this correctly for a one-sided power spectral density.

"Heinrich Acker" wrote in message <jg950d$g5k$1@newscl01ah.mathworks.com>...
> Hello Matlab users,
>
> I have used the code example in the documentation for the fft function
>
> http://www.mathworks.de/help/techdoc/ref/fft.html
>
> for my own program. The code claims to show an amplitude spectrum. The signal components shown in the example indicate this, and it can be easily tested with various numerical examples. There seems to be one exeption, though: The first element of the spectrum, 'the DC component'. The code
>
> Y = fft(y,NFFT)/L;
> f = Fs/2*linspace(0,1,NFFT/2+1);
> plot(f,2*abs(Y(1:NFFT/2+1)))
>
> appears wrong to me with respect to this element, because it is the only one that must not be multiplied by 2. I have tested this by adding a DC component to test signals, which appears with twice the amplitude in the spectrum. I would like to know what the experts say, because it could also be a matter of the definition of an amplitude spectrum, etc., etc.
>
> Heinrich

You're correct Heinrich, but even more than that, the Nyquist frequency only occurs once as well so that should not be doubled either. spectrum.periodogram in the Signal Processing Toolbox, gets this right:

Compare xdft and psdest.Data

n = 0:127;
x = cos(pi/2*n)+randn(size(n));
psdest = psd(spectrum.periodogram,x,'Fs',1,'NFFT',length(x));
%%%%
xdft = fft(x);
xdft = 1/length(x).*abs(xdft(1:length(x)/2+1)).^2;
% Do not double DC or the Nyquist
xdft(2:end-1) = 2*xdft(2:end-1);


Compare xdft and psdest.Data

norm(xdft'-psdest.Data)
max(abs(xdft'-psdest.Data))

Thanks you, Wayne, for clarifying this.

Regarding the Nyquist frequency: Isn't it wise to omit this point from an output called 'amplitude spectrum' anyway? The amplitude at this frequency can not be delivered, instead, we get only one component of the complex amplitude, which can be misleading, depending on the application.