Matlab plot log magnitude of chirp

Question

I create a logarithmic chirp exactly as found on the matlab help page.

t = 0:0.001:10;      % 10 seconds @ 1kHz sample rate
fo = 10; f1 = 400;   % Start at 10Hz, go up to 400Hz
X = chirp(t,fo,10,f1,'logarithmic');
figure(2);
spectrogram(X,256,200,256,1000,'yaxis');

Then I bring it to the frequency domain with the following code which works on other applications for me.

fft_prep = fftshift(fft(X));
fft_mag = abs(fft_prep);
pos_fft = fft_mag(1:ceil(length(fft_mag)/2));
db_fft = 20*log10(pos_fft);
figure(1);
plot(db_fft);

And I was surprised to see the following graph appear to be exciting 1kHz-5kHz:

I am not as familiar with the chirp function in matlab and was wondering if anyone saw something obvious that I am missing. Any other pointers are welcome.

Answer 1

Nothing wrong with the chirp function...

You just need to plot your db_fft against frequency values, and not vector indexes =).

plot(linspace(fo,f1,length(db_fft)), db_fft);

I also tested calculating the FFT of your signal using my other FFT methods and they too indicate a range between 0 and 400 Hz.

UPDATE:

IMO, I find it visually easier not to plot in dB or power (periodogram). Here's an excellent example and my goto method of calculating the FFT of a time-domain signal: mathworks.se/help/matlab/ref/fft.html

RESPONSE:

After some thinking, I agree I'm incorrect with my answer above, but not for the reason you say. The x-axis in the frequency domain should NOT go to the actual length of the chirp (or half, or dubbel or anything like that). The x-axis in the frequency domain should go to half of the sampling rate of the signal (Fs/2) and its your obligation to ensure you are sampling your signal with a sampling frequency twice that of the maxium frequency you wish/hope to resolve.

In other words, its incorrect to assume your FFT is of same/twice/half the length of your time domain signal because we can choose ANY number of frequency bins to represent the FFT in, and best practice is a length = N^2 (power of 2) for quick computation. Think about, why do you even need to know the time values when you calculate the FFT? You dont! You only need the sampling frequency (which should be set to Fs = 1000 btw, not Fs = 0.001).

My answer above then is incorrect, it should be:

plot(linspace(0, Fs/2, length(db_fft)), db_fft)

Instead of Fs/2, you have written length(t)/(2*Tfinal). it's (almost) the same value as Fs/2 but its not the correct way going about it =).

Here is my goto FFT method (values not in dB).

function [X,f] = myfft(x,Fs,norm)
    % usage: [X, f] = myfft(x,Fs,norm);
    %        figure(); plot(f,X)
    % norm: 'true' normalizes max(amplitude(fft))=1, default=false.
    if nargin==2
        norm=false;
    end
    L = length(x); NFFT = 2^nextpow2(L);
    f = Fs/2*linspace(0,1,NFFT/2+1);
    %f =0:(Fs/NFFT):Fs/2;
    X = fft(x,NFFT)/L; X = 2*abs(X(1:NFFT/2+1));
    if norm==true; X = X/max(abs(X)); end
end

And here's the resulting graph from [Xfft, f] = myfft(X,Fs); plot(f,Xfft); Notice that the return frequency bin vector has max(f) = Fs/2 in accordance to NyQuist theorem (any higher frequencies than Fs/2 cannot be resolved).