How to plot fourier series in matlab

Question

For homework (with no prior experience in matlab, guh.) I'm asked to do the following:

Plot the (2N + 1)-term approximation $$\sum\limits_{k=-N}^N{a_ke^{jk\omega_0t}}$$

where $a_k = \frac{\sin(k\omega_0T_1)}{k\pi}, k \neq 0 $, for N = 5, 21, 101, 501. Explain the behavior and verify the Gibbs phenomenon ($x(t)$ is a square wave). Assume $T = 16$ and $T_1 = 4$ and show the plots for one period from $\frac{-T}{2}$ to $\frac{T}{2}$

So now that that's out of the way ... my problem! I've reduced $a_k$ to $a_k = \dfrac{\sin(\frac{k\pi}{2})}{k\pi}, k \neq 0 $, since $\omega_0 = \frac{2\pi}{T}$ and I'm given $T$ and $T_1$.

So, substituting everything in, I get$$\sum\limits_{k=-N}^N{\dfrac{\sin(\frac{k\pi}{2})}{k\pi}e^{jk\frac{pi}{8}t}}$$

My issue is: how do I put that into Matlab and have it plot over $t$ and show the Gibbs phenomenon?

What I tried was

plot(t,(sum(sin(k*2*pi*0.25)./(k*pi))*exp(i*k*pi*t/8)), k = -5..5))

I just found symbolic variables, so I'll use those (I'm pretty sure I simplified wrong), but that doesn't resolve my issue of not knowing how to do the task in the first place.

I get an error re: k, "Error: The expression to the left of the equals sign is not a valid target for an assignment."

So how do I plot a Fourier series in Matlab? Conceptually it should be ... a series of 2k+1 vectors in t, summed together, but I don't understand how to do that, or how to plot the result, let alone use that to demonstrate the Gibbs phenomenon...

Answer 1

T = 16;
T1 = 4;
Nvec = [ 5, 21, 101, 501];

% The time samples at which we evaluate
% the signal.
dt = 1/100;
t = ( -T/2 : dt : +T/2 )';
Nt = length(t);

for N = Nvec

  k = ( -N : 1 : +N );

  % Define the sinc function that makes up the
  % series' weights. We need to catch the case
  % where the denominator is 0 and set it
  % to 1 (L'Hopital's rule).
  zinds = ( k == 0 );
  a = sin( k * 2*pi/T * T1 ) ./ (k*pi);
  a(zinds) = 1;

  % Define the terms in the approximation and sum
  % them. The below will make the kth term in the
  % sum the kth column in a matrix. Then we sum
  % the matrix across the columns.
  x = sum( repmat(a, Nt, 1)...
    .* exp(1j .* repmat(k,Nt,1) .* 2*pi/T .* repmat(t,1,length(k))), 2 );

  % Plot the resulting approximation.
  figure;
  plot( t, abs( x ), 'k' );
  title( sprintf( 'N = %d', N ) );

end