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).

Answer 2

I had a couple of errors, but not all of them were fixed. After trying the code out a little more, here is the solution I have come up with.

I added a few more variables to make the setup more modular.

Tfinal = 10;
Fs = 1/1000;
t = 0:Fs:Tfinal;      % 10 seconds @ 1kHz sample rate
fo = 10; f1 = 400;   % Start at 10Hz, go up to 400Hz
X = chirp(t,fo,Tfinal,f1,'linear');

When I plot the magnitude versus the linspace, the linspace needs to match up to the length of the actual chirp signal and not just from the low frequency to the high frequency. Because the vector t is of length 1000 and after the FFT we use the positive half, the chirp signal is 500 long and not 400, which the linspace up to f1 would suggest.

fft_prep = fftshift(fft(X));
fft_mag = abs(fft_prep);
pos_fft = fft_mag(ceil(length(fft_mag)/2)+1:length(fft_mag));
db_fft=20*log10(pos_fft);
figure(1);
plot(linspace(0,length(t)/(2*Tfinal),length(db_fft)), db_fft);

I also made the fix previously mentioned to get the positive side of the FFT instead of the negative. This plots:

This accurately graphs the chirp exciting 10Hz-400Hz. You can see it clearer by going to extreme cases and making it linear. Try a sampling frequency of 100, and a range from 10-25, without the linspace correction:

After the correction:

Answer 3

Incidentally, it may be useful to extract the original chirp amplitude from the FFT. For example if you are doing an audio sweep in frequency of some device, and you would like to know the amplitude and phase of the response at each frequency (like Pspice givs you). In short, the amplitude is a^2 = abs(FFT)^2 *4 *bandwidth of chirp /(Fs * N) where Fs is sample frequency and N is the number of points in the FFT. e.g. the bandwidth of a chirp from 200 to 400Hz is 200Hz.

If you want to know the dervition, start with Parseval's theorum: mean sqr of time signal = area under PSD. Thus, a^2/2 = sum(abs(FFT)^2) / N^2 where a is the peak amplitude of the swept frequency signal (e.g. chirp). If the chirp is linear rather than logarithmic in frequency distribution, then the FFT is flat, as in the above plots. THus we may replace the sum with Nb * abs(FFT)^2 / N^2 where Nb is the number of frequency bins that the chirp occupies, and abs(FFT) is the magnitude of the FFT, which is the same for all the frequency bins that the chirp occupies. Using the fact that the bandwidth of one frequency bin is Fs/N, we get that the bandwidth of the chirp is Nb * Fs /N. The result above is now easily derivable from here.

Answer 4

MATLAB's documentation about fft actually provides simple instructions that work generically for any choice of chirp (e.g., quadratic or linear):

Fs = 1000;            % Sampling frequency                    
T = 1/Fs;             % Sampling period       
L = 1500;             % Length of signal
t = (0:L-1)*T;        % Time vector

Example of signal to be analysed:

S = 0.7*sin(2*pi*50*t) + sin(2*pi*120*t);
X = S + 2*randn(size(t));

Calculate FFT:

Y   = fft(X);       % Calculate FFT
P2  = abs(Y/L);
P1  = P2(1:L/2+1);
P1(2:end-1) = 2*P1(2:end-1);

The frequency vector is calculated as

f = Fs*(0:(L/2))/L;

If you then want to generate a Bode magnitude plot, you first should convert frequency to [rad/s] and the outcome of the fft to [dB]:

fRad = f*2*pi;
Pdb = 20*log10(P1);

and then make a Bode plot (I suggest to use scatter, given the potential noisy nature of the outcome)

figure
scatter(fRad,Pdb)
set(gca,'xscale','log') 
grid on; grid minor
xlabel('Frequency [rad/s]')
ylabel('Magnitude [dB]')