Search code examples
matlabsignal-processingfftdft

There are different results when I take FFT sinc(t)


The first situation:

%% Analytical calculation
syms t f
h = @(f) int(sinc(t)*exp(-1i*2*pi*f*t),t,-10,10);
subplot(2,1,1); fplot(real(h(f)),[-3 3]); grid; ylim([-3 3]);
subplot(2,1,2); fplot(imag(h(f)),[-3 3]); grid; ylim([-3 3]);

The second situation:

%% Numerical calculation
N=100;
T = 10;
t = (2*(0:N-1)/N-1)*T;

x = sinc(t);
y = (fftshift(fft(x)));
figure;
subplot(2,1,1); plot((real(y)));
subplot(2,1,2); plot(imag(y));

Why is the result different? Why isn't rectangular impulse in the second case? What is the difference between abs and real and what should I use correctly?

Analytical calculation

Numerical calculation


Solution

  • Short answer

    In the second part, replace your line

    x = sinc(t);
    

    by

    x = fftshift(sinc(t));
    

    so that the sinc function is centered at the time origin, which is the first sample.

    Long answer

    The DFT (FFT) assumes that the time origin is at the first sample, regardless of the time axis t that you are define. Note that in your code

    t = (2*(0:N-1)/N-1)*T;
    x = sinc(t);
    y = (fftshift(fft(x)));
    

    the fft function knows nothing about your t variable. It simply considers the time axis to be [0, 1, ..., numel(x)-1]. Since your x is defined as a sinc that is centered within the observation window (its main lobe in the middle), the DFT interprets that as a sinc that is not centered around the origin, but rather it is (circularly) shifted by half the observation window.

    As you may know, a shift in time corresponds to multiplying by an exponential in the frequency domain. For the DFT, this property applies for circular shifts.

    Therefore, the plot you get differs from the one you expected in that it has been multiplied by an exponential at half the sample rate, which reduces to [1, -1, 1, -1, ...]. Note how the values in your plot are alternatively positive and negative.

    To solve this, you need to define the sinc centered at the origin, that is, at the first sample. So the peak and the right half of the main lobe will be in the leftmost region of the observation window, and the left half of the lobe will appear in the rightmost region. This can be done very simply by applying fftshift to your x:

    x = fftshift(sinc(t));
    plot(x)
    

    enter image description here

    Then the DFT is

    y = (fftshift(fft((x))));
    subplot(2,1,1); plot((real(y)));
    subplot(2,1,2); plot(imag(y));
    

    enter image description here