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

Fourier Series Expansion using Mathematica

  1. Nov 30, 2011 #1

    I recently learned about the Fourier Series and how it can be used decompose a periodic signal into a sum of sinusoids. I can calculate all the coefficients by hand, but I wanted Mathematica to do that for me. I attempted to write a code, and it does give the desired output.

    I wanted to animate the output, in the sense that Mathematica would display more and more corrections over time. So I want to start by displaying 1 partial sum, and over time display many partial sums. I looked into the Animate function to do this, but ran into some problems as it requires an index to iterate over. Naturally, I chose the number of partial sums included as the iterator. So instead of having the number of partial sums to include as an integer, I make it a variable. But doing so leaves Mathematica in some infinite loop. Below is my code. I tried to include as many comments as possible, and would appreciate any help.

    Code (Text):
    (*Define basis functions of sin and cosine*)
    cosBasis[m_, t_, T_] := Cos[2*Pi*m*t/T];
    sinBasis[m_, t_, T_] := Sin[2*Pi*m*t/T];

    (*Define inner product*)
    fourierIP[f_, g_] := Integrate[f*g, {t, -Pi, Pi}]

    (*Define function to calculate coeffecients of sin and cosine,
    depends on (f[t],period,m)*)
    fourierSinCoeff[func_, T_, m_] :=
      fourierIP[func[t], sinBasis[m, t, T]]/
       fourierIP[sinBasis[m, t, T], sinBasis[m, t, T]];
    fourierCosCoeff[func_, T_, m_] :=
      fourierIP[func[t], cosBasis[m, t, T]]/
       fourierIP[cosBasis[m, t, T], cosBasis[m, t, T]];

    (*Compute M^th partial sum of fourier coeffecients*)
    fourierSeries[func_, t_, T_, M_] :=
      Sum[fourierCosCoeff[func, T, m]*cosBasis[m, t, T], {m, 0, M}] +
       Sum[fourierSinCoeff[func, T, m]*sinBasis[m, t, T], {m, 1, M}];

    const = 0.5;(*height of square wave*)
    squareWave[t_] :=
      Piecewise[{{const, 0 < t <= Pi}, {0,
         Pi < t <= 2*Pi}}];(*construct square wave*)
    periodicExtension[func_, nPeriods_] :=
     Sum[func[t + 2*Pi*n], {n, -nPeriods,
       nPeriods}];(*extend square wave over multiple periods**)
    plot1 = Plot[periodicExtension[squareWave, 4], {t, -4*Pi, 4*Pi},
      PlotRange -> {{-4*Pi, 4*Pi}, {-1, 1}},
      ExclusionsStyle -> Dotted](*display square wave*)
    fourierCosCoeff[f, 2*Pi, m];
    fourierSinCoeff[f, 2*Pi, m];
    output = fourierSeries[squareWave, t, 2*Pi,
      9];(*include 9 partial sums*)
    plot2 = Plot[output, {t, -4*Pi, 4*Pi}];
    Show[plot1, plot2](*displays plot of square wave and fourier wave together*)
    (*The above code works just fine, the moment I switch the
    number of partial sums to include into a variable, I get problems*)

    (*Now to animate over many partial sums*)
    output1 = fourierSeries[squareWave, t, 2*Pi, m];
    Animate[output1, {m, 1, 3, 1}]
  2. jcsd
  3. Nov 30, 2011 #2
    Never mind, I got it.

    The problem seems to be that mathematica treats the iterator as a general variable (ie. one that can be numbers,strings,etc.) and so takes what seems like forever to compute the expression for the fourier series. By specifying the iterator as an integer, the process takes a few seconds.

    Code (Text):
    Compile[{m, _Integer}, m];
    output1 = fourierSeries[squareWave, t, 2*Pi, m];
    Animate[Plot[output1 /. m -> k, {t, -4*Pi, 4*Pi}], {k, 1, 100, 1}]
Share this great discussion with others via Reddit, Google+, Twitter, or Facebook