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]);
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')
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;
Y = fft(y)/Fs; % Normalize FFT output % to take into account that Fs is 5000, not 1
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:
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:
Just make sure you normalize the amplitudes in "Y" by Fs as I indicate above!
Also, to get an even time series change:
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.
You might want to check out