Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

FFT in Matlab

  1. Jul 31, 2011 #1
    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 [tex]cos(2*\pi*x*f_0)[/tex], which should be [tex]\frac{\delta (f + f_0) + \delta (f - f_0)}{2}[/tex]

    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 [tex]H \left( \frac{n}{NT} \right) = \sum_{k=0}^{N-1} h(kT) e^{-i2\pi n k / N} [/tex] for n = 0, ..., N - 1.
     
  2. jcsd
  3. Jul 31, 2011 #2

    Pythagorean

    User Avatar
    Gold Member

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

     
     
  4. Jul 31, 2011 #3
    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.
     
  5. Jul 31, 2011 #4

    Pythagorean

    User Avatar
    Gold Member

    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.
     
  6. Jul 31, 2011 #5
    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....
     
  7. Aug 1, 2011 #6

    Pythagorean

    User Avatar
    Gold Member

    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]
     
  8. Aug 1, 2011 #7
    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
  9. Aug 2, 2011 #8

    Pythagorean

    User Avatar
    Gold Member

    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.
     
  10. Aug 2, 2011 #9
    Not sure what you mean by internally. My text, matlab, and wikipedia all use DFT convention [tex]H_n = \sum_{k=0}^{N-1} h(kT) e^{-i 2 \pi n k / N}[/tex]. Where T is sampling interval and N is number of points sampled.
     
  11. Aug 3, 2011 #10

    Pythagorean

    User Avatar
    Gold Member

    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.
     
  12. Aug 3, 2011 #11
    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).
     
  13. Aug 3, 2011 #12
    See attached picture.

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

    Attached Files:

  14. Aug 3, 2011 #13

    Pythagorean

    User Avatar
    Gold Member

    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
     
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook




Similar Discussions: FFT in Matlab
  1. Fft in matlab (Replies: 0)

  2. Matlab (fft) (Replies: 1)

  3. Matlab fft help (Replies: 10)

  4. MATLAB FFT question (Replies: 2)

Loading...