# FFT in Matlab

1. Jul 31, 2011

### 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.

Code (Text):

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.

2. Jul 31, 2011

### 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:

Code (Text):

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));

3. Jul 31, 2011

### 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.

4. Jul 31, 2011

### 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.

5. Jul 31, 2011

### 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....

6. Aug 1, 2011

### 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.)

Code (Text):

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]

7. Aug 1, 2011

### PhDorBust

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

Last edited: Aug 1, 2011
8. Aug 2, 2011

### 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.

9. Aug 2, 2011

### PhDorBust

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.

10. Aug 3, 2011

### Pythagorean

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:

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.

11. Aug 3, 2011

### PhDorBust

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).

12. Aug 3, 2011

### PhDorBust

See attached picture.

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

#### Attached Files:

• ###### question.png
File size:
7.7 KB
Views:
135
13. Aug 3, 2011

### Pythagorean

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