Question about fft

I was doing a small program to graph the fft of a discrete signal, and to validate it started trying some functions whose Fourier transforms are known, I initially tried "sin(2*pi*w*t)" and "cos(2*pi*w*t)" and the result was, as expected, a "spike" at -w and w (when plotting the absolute value of the fft), then I tried the rectangle function (I think it's also known as the top hat function), and i expected to get something similar to "sinc(x)" as an answer, but instead i get the following graph: Not_Sinc, the original signal is Rect. If I zoom the image around 0 I can see something that looks like sinc, but there seems to be a delta in 0 that should not be there.

I tried increasing the sampling frequency and the results get worse.

Thanks in advance for your answers!

The code I used is here:

clear all
close all
clc

Fs=5000;   %[Hz] Sampling frequency
x=(-1:1/Fs:1)';  %Time
y=zeros(size(x));  %Preallocate vector
for i=1:length(x)  %Awfully ineficcient way to create rect function
    if x(i)>=-0.5 && x(i)<=.5
        y(i)=1;
    end
end

%Plot of the signal
plot(x,y)
xlabel('time(s)');
ylabel('Amplitude');
axis([-Inf Inf -0.1 1.1]);

Y=fft(y);
Y=fftshift(Y);

d=Fs/length(Y(:,1));  %delta for the frequency vector
F=-Fs/2:d:Fs/2-d;     %limits due to Nyquist

figure
plot(F,(sqrt(Y.*conj(Y))))
xlabel('Frequency [Hz]')
ylabel('Power')

Answer 1

[Dr. Seis]

The way it is implemented above, the value at 0 frequency is just the sum of your signal from -1 seconds to 1 second (you have about 5000 samples equal to 1 and about 5000 samples equal to 0)... so the expected amplitude should be about 5000.

A couple of suggestions that won't change much:

y = zeros(size(x));
y(x(i) >= -0.5 & x(i) <= 0.5) = 1;

and

Y = fft(y)/Fs; % Normalize FFT output 
               % to take into account that Fs is 5000, not 1

and

abs(Y) is equivalent to sqrt(Y.*conj(Y))

EDIT May 09, 2012

To plot something that looks more like a sinc, change your line:

plot(F,(sqrt(Y.*conj(Y))))

to

plot(F,real(Y))

and narrow your x-limits to +/- 6 (i.e., "xlim([-6 6])")

To confirm it is what you expect, create the sinc-like function you expect. From the "DFT pairs" link below, this would look something sort of like:

Z = zeros(size(Y));
Z(5001)=1;
Z(5002:5051) = -1*sin(pi*10000*(1:50)/20001)./(10000*sin(pi*(1:50)/20001));

Not sure why I needed the -1, but the positive frequency amplitudes seem to line up pretty well with what the expected shape is when you plot:

plot(F,real(Y),F,Z);

Just make sure you normalize the amplitudes in "Y" by Fs as I indicate above!

Also, to get an even time series change:

x=(-1:1/Fs:1)';

to

x=(-1:1/Fs:1-1/Fs)';

Answer 2

[Esteban Zegpi]

I know the way to create the vector isn't very elegant, I didn't care because I will read the input from files later on, so it was a non issue. And about the normalization, that will not change the fact that the value at f=0 should be the area of the square and not the amount of samples, as you can see here DFT pairs, the "discrete sinc" should be the Discrete Fourier Pair of the rect function.

Regards

[Dr. Seis]

In fact... it does change the fact.

What is the area of the square? Answer: 1

What is the value of the FFT at zero frequency when you normalize the way I told you? Answer: 1

The Fourier transform is the integral of your data in the time domain multiplied with some exponential function, yes? When you evaluate the forward Fourier transform at 0 frequency, the exponential function disappears (i.e., exp(0) = 1). Since Matlab assumes that Fs = 1 when it does the FFT, the Fourier transform of your time domain data (i.e., "y") at zero frequency ends up being equivalent to "sum(y)". In order to get an accurate representation of the amplitude in the frequency domain, you must normalize by the true Fs to account for the fact that width between samples (area under curve is length x width) is not Fs=1. Otherwise, as in your case, the area will be off by a factor of exactly Fs.

[Esteban Zegpi]

What I ment to say was that normalizing by the sampling frequency would not change the fact that the "sinc" function would still have a spike in f=0, and look nothing like a sinc. I realize that the normalization you proposed is correct.

[Dr. Seis]

It doesn't look like a sinc function because you are effectively plotting the "abs(Y)" (no negative amplitudes)... try plotting "real(Y)". It doesn't look exactly the same (but close) as the example here:

http://en.wikipedia.org/wiki/Fourier_transform

Scroll down about 1/5 from the top is an example of a Top-hat function and its' Fourier Transform. Also stick an "xlim([-6 6])" below where you plot the frequency domain stuff.

Answer 3

[Daniel Shub]

You might want to check out

http://blogs.mathworks.com/steve/category/fourier-transforms/

[Esteban Zegpi]

Thank you, I'll check this out!