CODEWITHSUNDEEP
matlabfft
LudovicG
LudovicG

Reputation: 25

Trying to plot the fft of a sinc function

I am trying to plot the fft of a set of data I have. These data form a nearly perfect sinc function. Here is the data of which I am trying to plot the fft: Here is the data of which I am trying to plot the fft.

I know the fft of a sinc function should look like kind of a step function. However, the results I get are nowhere near that. Finding the fft in itself is super easy, but I think my mistake is when I try to compute the frequency axis. I have found several methods online, but so far I have not been able to make it work. Here is my code:

sampleRate = (max(xdata) - min(xdata))/length(xdata);
sampleN = length(xdata);
y = fft(ydata, sampleN);
Y = y.*conj(y)/sampleN;
freq = (0:1:length(Y)-1)*sampleRate/sampleNumber;
plot(freq, Y)

I have found pretty much all of that online and I understand pretty much none of it (which might be why it's not working...)

Zoom on what I get using that code: Zoom on what I get using that code

It now seems to be working! This is what I get when I subtract the mean: This is what I get when I subtract the mean

Upvotes: 1

Views: 1285

Answers (1)

Cris Luengo
Cris Luengo

Reputation: 60554

What you see here is the zero frequency being much, much larger than everything else. Plot with plot(freq,Y,'o-') to prove that the shape you see is just the linear interpolation between two samples.

The zero frequency is the sum of all samples. Because the mean of your signal is quite a bit larger than the amplitude, the zero frequency dwarfs everything else. And because you are plotting the power (absolute square of the DFT), this difference is enhanced even more.

There are two simple solutions:

  1. Plot using logarithmic y-axis:

    plot(freq, Y)
set(gca,'yscale','log')

  2. Subtract the mean from your signal, remove the zero frequency, or scale the y-axis (these are all more or less equivalent):

    y = fft(ydata-mean(ydata), sampleN);

    or

    y(1) = 0;

    or

    plot(freq, Y)
set(gca,'ylim',[0,max(Y(2:end))]);

Upvotes: 1

Related Questions