Getting Close with Matlab FFT for Fourier Transform Table Reconstruction

[PhDorBust]

I'm attempting to use Matlab fft functionality to reconstruct Fourier transform tables in my textbook (brigham), but to little success.

Here is code to take the Fourier transform of cos(2*\pi*x*f_0), which should be \frac{\delta (f + f_0) + \delta (f - f_0)}{2}

I can *almost* get it, but not quite. See spikes at +/- .001, when should be at +/- .1 Any help would be appreciated.

 x = 0:.1:9.9;
%N=100, T=10
 y = cos(x*2*pi/10);
 Y = fft(y);

X = (0:99)/1000 - 50/1000;
Y = fftshift(Y);
plot(X,real(Y));

Using H \left( \frac{n}{NT} \right) = \sum_{k=0}^{N-1} h(kT) e^{-i2\pi n k / N} for n = 0, ..., N - 1.

 

[Pythagorean]

you defined X arbitrarily. Your frequency resolution is set by the length of your original signal (x), and your sample rate (fs), which is of course, 1/dx (dx being the "period" of each sample).

Your frequency range is set by the nyquist frequency, which is half the sampling rate.

So try this:

dx = .1;
 x = 0:dx:9.9;

fs = 1/dx;

 y = cos(x*2*pi/10);
 Y = fft(y);

X = -fs/2 : fs/length(x)  : fs/2 - fs/length(x)
Y = fftshift(Y);
plot(X,real(Y));

 

[PhDorBust]

You're correct, I wasn't setting the frequency domain properly. I.e. T is .1, not 10.. Thank you.

Also, that is nicer notation.

 

[Pythagorean]

Ah, yes, I see your range vector was actually quite intended for something.

A good sci/eng coding philosophy is to try to not hard-code numbers; try to use lots of intuitive, explanatory variable names so you can keep track of what's going on.

 

[PhDorBust]

Question 2: The amplitudes of the impulses should be 1/2.

For band-limited, periodic functions such as cosine, the exact Fourier transform should be given by multiplying the DFT by T (the sampling interval). But I always get 5 for amplitudes, when I want .5...

 

[Pythagorean]

yeah, we should normalize by the length of the signal from which the fft is taken (which should be the length of fft too, as long as you don't designate that. It's best not too, you can pad with zeros yourself (and always should to make your vectors of length 2^n for efficiency.)

dx = .1;
 x = 0:dx:9.9;

fs = 1/dx;
F1 = .1;
F2 = .2;
 y = cos(x*2*pi*F1);
 Y = fft(y);

X = -fs/2 : fs/length(y)  : fs/2 - fs/length(y)
Y = fftshift(Y);
[B]plot(X,real(Y)/length(y));[/B]

 

[PhDorBust]

yeah, we should normalize by the length of the signal from which the fft is taken (which should be the length of fft too, as long as you don't designate that. It's best not too, you can pad with zeros yourself (and always should to make your vectors of length 2^n for efficiency.)

dx = .1;
 x = 0:dx:9.9;

fs = 1/dx;
F1 = .1;
F2 = .2;
 y = cos(x*2*pi*F1);
 Y = fft(y);

X = -fs/2 : fs/length(y)  : fs/2 - fs/length(y)
Y = fftshift(Y);
[B]plot(X,real(Y)/length(y));[/B]

The way you are doing it appears to be correct. But my textbook (brigham) says "If we desire to compute the Fourier transform by means of the discrete Fourier transform, it is necessary to multiply the discrete time function by the factor T", where T is the sampling interval (dx in your code). This is true for band-limited periodic functions.

See http://en.wikipedia.org/wiki/Relati...eries#DFT_versus_continuous_Fourier_transform

 

[Pythagorean]

depending how you do the transform internally, you have to normalize it externally. There's not really a single convention for how to do this.

Your book is probably correct in their treatment, but that's a different treatment than matlab's.

 

[PhDorBust]

depending how you do the transform internally, you have to normalize it externally. There's not really a single convention for how to do this.

Your book is probably correct in their treatment, but that's a different treatment than matlab's.

Not sure what you mean by internally. My text, matlab, and wikipedia all use DFT convention H_n = \sum_{k=0}^{N-1} h(kT) e^{-i 2 \pi n k / N}. Where T is sampling interval and N is number of points sampled.

 

[Pythagorean]

Not sure what you mean by internally.

If you're just starting programming, then you probably can't appreciate that going form theory to code isn't always as straightforward as you'd think. There's lots of different ways to do things and a lot of different things that customers want out of a function/program.

If you want to be sure how your fft works, write one yourself. Don't use matlab's.

If you want a more in-depth discussion on fft normalization conventions:

http://www.mathworks.com/matlabcentral/newsreader/view_thread/152029

Remember that the fft is NOT inherently physical. YOU have to implement the physics. It can be normalized differently depending on the question you're asking. But in order to do that, you have to know what the results mean when they come out of the function that you're using so that you can normalize to your physical system.

 

[PhDorBust]

If you're just starting programming, then you probably can't appreciate that going form theory to code isn't always as straightforward as you'd think. There's lots of different ways to do things and a lot of different things that customers want out of a function/program.

If you want to be sure how your fft works, write one yourself. Don't use matlab's.

If you want a more in-depth discussion on fft normalization conventions:

http://www.mathworks.com/matlabcentral/newsreader/view_thread/152029

Remember that the fft is NOT inherently physical. YOU have to implement the physics. It can be normalized differently depending on the question you're asking. But in order to do that, you have to know what the results mean when they come out of the function that you're using so that you can normalize to your physical system.

If you say so, but all the Cooley-Tukey in FFT does is evaluate the sum I wrote out in a speedy way. I'm not really normalizing... just want to match the DFT result to FT which should be scaling by T (see wiki link I posted).

 

[PhDorBust]

See attached picture.

When taking the DFT, why do we have to take the time waveform to be (b) rather than (a)?

 

[Pythagorean]

just want to match the DFT result to FT which should be scaling by T (see wiki link I posted).

Your wiki link says that DFT = (1/T)*FT. You want DFT to look like FT, so you solve for FT:

FT = T*DFT

For a discrete n sec. signal, the sampling rate T is inversely proportional to the N, the number of samples. (n = T*N)

so T ~ 1/N