Sine and Square Waves for St. Patrick's Day

This notebook is part of a guest lecture I gave in our introductory signal processing class, ELEG 305, on St. Patrick's Day.

The command below loads numpy (the vector processing library) and matplotlib (Matlab-style plotting) and tells the notebook to draw the plots in place. Many Ipython notebooks start with this command.

In [53]: 
%pylab inline

Populating the interactive namespace from numpy and matplotlib


WARNING: pylab import has clobbered these variables: ['f']
`%matplotlib` prevents importing * from pylab and numpy

Sine Wave

Let's plot a sine wave. The Nyquist theorem says we have to sample at least twice per period. Let's try it.

In [54]: 
f = 1
Fs = 2.1 #more than twice per period
t = linspace(0,3,3*Fs+1)
y = sin(2*pi*f*t)
plot(t,y)
Out[54]:
[<matplotlib.lines.Line2D at 0x10c7a6410>]

That didn't work. Clearly the plotting routine does not understand the Nyquist criterion (neither, by the way, does Matlab's plotting routine).

We need to increase the sampling rate.

In [55]: 
f = 1
Fs = 100 #a lot more than twice per period
t = linspace(0,3,3*Fs+1)
y = sin(2*pi*f*t)
plot(t,y)
Out[55]:
[<matplotlib.lines.Line2D at 0x10d9765d0>]

That's better. They look like sinusoids.

Lissajous figures

Various plots of sines and cosines on both axes.

In [56]: 
t = linspace(0,1,201)
wt = 2*pi*t
plot(cos(wt), sin(wt))
axes().set_aspect('equal')

In the plots below, play with the frequencies and the phase shift to change the picture. There are many interesting Lissajous figures.

In [57]: 
plot(cos(4*wt),cos(3*wt+pi/5))
axes().set_aspect('equal')

Let's improve the picture a bit by increasing the box size a little.

In [58]: 
X=cos(2*wt)
Y=cos(3*wt+pi/4)
axes().set_aspect('equal')
xlim(X.min() * 1.1, X.max() * 1.1)
ylim(Y.min() * 1.1, Y.max() * 1.1)
plot(X,Y)
Out[58]:
[<matplotlib.lines.Line2D at 0x10db98350>]

With a small change, we can generate other shapes. Here we want successive time values to be relatively far apart so the points on the circle jump around.

Let's also eliminate the box.

In [70]: 
tt=0.4*arange(0,6)
plot(cos(2*pi*tt),cos(2*pi*tt+pi/2))
axes().set_aspect('equal')
axes().axis('off') #comment out to get box back
Out[70]:
(-1.0, 1.0, -1.0, 1.0)

Square Waves

The Fourier series for a square wave is the following:

\[sq(t) = \frac{4}{\pi}\sum_{\text{$k=1$, $k$ odd}}^\infty \frac{\sin(2\pi k f t)}{k} = \frac{4}{\pi} \left( \sin(2\pi f t) + (1/3) \sin(2\pi 3 f t) + (1/5)\sin(2\pi 5 ft) + \ldots\right)\]

The various frequencies are called harmonics. The first harmonic is also called the fundamental. Notice this series has only odd terms. I.e., all the even harmonics are 0. (Other Fourier series may have non-zero even terms.)

Let's look at a some approximations to the square wave.

In [60]: 
t = linspace(0,2,201) #two cycles, 100 points each
y1 = (4/pi)*sin(2*pi*f*t)
y3 = y1 + (4/(3*pi))*sin(2*pi*3*f*t)
y5 = y3 + (4/(5*pi))*sin(2*pi*5*f*t)
sqwave = sign(sin(2*pi*f*t)) #an actual square wave
plot(t,y1, t,y3, t,y5, t, sqwave)
Out[60]:
[<matplotlib.lines.Line2D at 0x10c7b4410>,
 <matplotlib.lines.Line2D at 0x10c7b49d0>,
 <matplotlib.lines.Line2D at 0x10c714450>,
 <matplotlib.lines.Line2D at 0x10c7c0550>]

We can clearly see we are approaching a square wave. Let's add a bunch more terms, but let's also do it more efficiently.

In [61]: 
sq = zeros(len(t)) #preallocate the output array
for h in arange(1,25,2):
    sq += (4/(pi*h))*sin(2*pi*f*h*t)

plot(t,sq, t,sqwave)
Out[61]:
[<matplotlib.lines.Line2D at 0x10db36750>,
 <matplotlib.lines.Line2D at 0x10db369d0>]

Those little peaks near the discontinuities are called Gibbs Phenomena. As long as we have a finite number of terms, no matter how many, we will have Gibbs phenomena.

Note, if we change the sine to cosine the final picture looks a lot different. Phase matters.

In [62]: 
not_sq = zeros(len(t))
for h in arange(1,15,2):
    not_sq += (4/(pi*h))*cos(2*pi*f*h*t)

plot(t,not_sq, t,sqwave)
Out[62]:
[<matplotlib.lines.Line2D at 0x10d9c7090>,
 <matplotlib.lines.Line2D at 0x10d9c7310>]

C. Boncelet, 17 March 2014, rev. 30 April 2014

In [62]: