Fun with Sinusoids -- And a Little Digital Music¶
Talk given at the Toyota-UD Applied Math Initiative, 24 July 2014.
Charles Boncelet
Professor, Electrical & Computer Engineering
University of Delaware
Newark DE 19711
boncelet@udel.edu
302-831-8008Updated 5 August 2014.
About this presentation -- Python and IPython¶
This presentation is an IPython notebook. The various calculations, graphs, and audio outputs are created using Python and some of its libraries (principally numpy and matplotlib).
The presentation is available in two forms:
An IPython notebook file with extension .ipynb. This is an executable file. It can be run in an IPython server.
An HTML file. This can be displayed in any browser, but it is frozen and cannot be edited.
Python is an easy to learn, but still very powerful language. It is the most popular first language in computer science programs at universities in the USA. At Delaware, engineering students use Python in their first programming course.
This notebook serves as an introduction to numerical computing with Python. Think of Python as a "super-calculator".
I recommend you install one of the following (free) packages:
Canopy. This is the package I use. It comes with lots of libraries (more than 100). It's free to academic users, but you have to ask for an academic license from Enthought.
Anaconda. This is newer than Canopy (I started with Canopy before Anaconda was available). Lots of people recommend Anaconda. It is free without any license issues and has lots of libraries. (If starting today, I'd try Anaconda.)
Unfortunately it does not seem possible (at least on my Mac) to run both Canopy and Anaconda on the same machine (to compare them). They seem to conflict with each other.
For this talk, I had to upgrade the IPython library to the 3.0.dev version. The current release, version 2.1, does not have the Audio command we use to run audio snippets in the browser. Hopefully, the 3.0 version is released soon and quickly incorporated into Canopy and Anaconda.
Miscellaneous Ipython and Python commands¶
%pylab inline
from IPython.display import Audio, Image
from __future__ import division #so 1/2 = 0.5, not 1/2=0
#ignore the error below about clobbering 'f', unimportant
Sinusoids¶
Sinusoids are used for all sorts of signal models in Electrical & Computer Engineering (ECE). It is hard to overstate their importance in ECE.
The basic sinusoid is a sine or cosine function over time.
\[ X(t) = A \cos(\omega t+\theta)\]
where
- \(A\) is the amplitude of the sinusoid.
- \(\omega\) is the frequency (in radians per second). \(\omega=2\pi f\) where \(f\) is the frequency in cycles per second (Hertz, abbreviated Hz).
- \(f=1/T\) where \(T\) is the period of the sinusoid, i.e., the sinusoid repeats every \(T\) seconds.
- \(\theta\) is a phase offset, in radians.
Note,
\[\cos(\omega t+\theta) = \sin(\omega t +\theta+\pi/2)\]
So we can use either \(\cos\) or \(\sin\) as our basic function.
Image('sinusoid-pic.png')
Image('Liver-sinusoid_colour_cartoon.jpg')
f=1
w=2*pi*f
t=linspace(0,3,25)
scatter(t,cos(w*t))
scatter(t,cos(w*t)) #dots
plot(t,cos(w*t)) #connect the dots
t=linspace(0,3,151) #lots of points
scatter(t,cos(w*t))
plot(t,cos(w*t))
plot(t,cos(w*t),'b',label='cos')
plot(t,sin(w*t),'r',label='sin')
legend()
Lissajous Figures¶
Interesting plots of one sinusoid on the X axis and a different one on the Y axis. Vary the frequencies and phase to get a variety of interesting pictures.
wt = linspace(0,2*pi,300)
X = cos(wt)
Y = sin(wt)
xlim(X.min() * 1.1, X.max() * 1.1)
ylim(Y.min() * 1.1, Y.max() * 1.1)
plot(X,Y)
axes().set_aspect('equal') #gives circle, not ellipse
X=cos(2*wt)
Y=cos(3*wt+pi/4)
xlim(X.min() * 1.1, X.max() * 1.1)
ylim(Y.min() * 1.1, Y.max() * 1.1)
plot(X,Y)
axes().axis('off') #eliminate the box
axes().set_aspect('equal')
Frequencies High and Low¶
Frequencies in applications vary across many orders of magnitude:
- Audio:
- Audible frequencies from 20 Hz to 20,000 Hz
- A above middle C: 440 Hz
- Ultrasound Imaging: about 2 MHz
- Power:
- 60 Hz in USA (e.g., hum in audio)
- 50 Hz in many other countries
- Radio:
- AM Broadcast: about 1 MHz
- FM Broadcast: about 100 MHz
- Cell Phones: about 1 GHz
- Microwave Oven: 2.45 GHz
- UD Research: up to 100 Ghz
Complex Numbers and Euler's Formula¶
In analyzing sinusoids, we make great use of complex numbers.
\[z = x + jy = |z|\angle{\theta}\]
where \(|z| = \sqrt{x^2+y^2}\), \(\theta = \tan^{-1}(y/x)\), \(x=\operatorname{Re}(z)\), and \(y=\operatorname{Im}(z)\).
\[\begin{align} z_1 + z_2 &= (x_1 + x_2) + j(y_1+y_2)\\ z_1 z_2 &= (x_1x_2 -y_1y_2) + j(x_1y_2 + x_2y_1)\\ &= |z_1||z_2| \angle(\theta_1+\theta_2) \end{align}\]
We use this simple formula below to do a hard computation:
\[\operatorname{Re}(z_1) + \operatorname{Re}(z_2) = \operatorname{Re}(z_1+z_2)\]
Perhaps the most important single formula in all of mathematics is Euler's Formula:
\[e^{j\theta} = \cos(\theta) + j\sin(\theta)\]
Using Euler's formula, we can also write a complex number as
\[z = |z| e^{j\theta}\] where \(\theta=\angle z\) is the angle of \(z\).
With \(\theta=\pi\), Euler's formula can be rearranged as
\[ e^{j\pi}+1=0\]
Thus combining the five most important mathematical values into one simple equation.
Use Complex Algebra -- Not Trigonometry!¶
While we use sinusoids in lots of areas, we don't use trigonometry. We use complex variables and Euler's formula.
Euler's formula gives us simple formula for sinusiods:
\[e^{j\omega t} = \cos(\omega t) + j\sin(\omega t)\]
Using Euler's formula, we can write sinusoids two different ways:
\[\cos(\omega t) = \operatorname{Re}(e^{j\omega t}) = \frac{e^{j\omega t} + e^{-j\omega t}}{2} \]
Using the \(\operatorname{Re}(e^{j\omega t})\) is especially handy when we add sinusoids together and the second formula when multiplying sinusoids.
Here's an example of adding two sinusoids (with the same frequency):
\[\begin{align} \cos(\omega t) + 2\cos(\omega t + \pi/4) &= \operatorname{Re}(e^{j\omega t}) + \operatorname{Re}(2 e^{j\omega t+j\pi/4})\\ &= \operatorname{Re}(e^{j\omega t}+2 e^{j\omega t+j\pi/4})\\ &= \operatorname{Re}((1+2 e^{j\pi/4})e^{j \omega t}) \end{align}\]
Now add the two complex numbers together and convert to polar notation:
z=1+2*exp(1j*pi/4) #use '1j' for the imaginary number
z
abs(z),angle(z)
Thus (with rounding), \[\begin{align}\cos(\omega t) + 2\cos(\omega t + \pi/4) &= \operatorname{Re}(2.80e^{j0.53}e^{j\omega t})\\ &= 2.80\operatorname{Re}(e^{j\omega t+j0.53})\\ &= 2.80 \cos(\omega t+0.53)\end{align}\]
Adding two sinusoids with the same frequency results in a single sinusoid of the same frequency.
Beat Frequency¶
Add two sinusoids with different frequencies:
\[\cos(\omega_1 t) + \cos(\omega_2 t) = 2 \cos\big((\omega_1-\omega_2)t/2\big)\cos\big((\omega_1+\omega_2)t/2\big)\]
\(\omega_1-\omega_2\) is the beat frequency.
Guitar tuners adjust one note until the beat frequency against a reference is 0.
Fs = 20050
f1 = 440
f2 = 441
t = linspace(0, 4, 4*Fs)
sound = cos(2*pi*f1*t) + cos(2*pi*f2*t)
Audio(sound, rate=Fs)
Frequency Shifting¶
Multiplying two sinusoids results in a frequency shift:
\[\cos(\omega_1 t) \cos(\omega_2 t) = \frac{\cos((\omega_1+\omega_2)t)}{2} + \frac{\cos((\omega_2-\omega_1)t)}{2}\]
This is how radio transmits a signal. \(\omega_1\) represents the sound (typically a small frequency) and \(\omega_2\) the radio carrier frequency (typically much higher).
Music Theory¶
Music frequencies are exponential. Each octave is a doubling in frequency.
E.g., Note A can be 55 Hz, 110 Hz, 220 Hz, 440 Hz, 880 Hz, etc.
\[A = 440*2^{m-4}\]
where \(m\) is octave number. \(m=4\) is the standard middle octave.
Music Scales¶
The conventional scale in western music consists of semitones \(2^{1/12}\) apart.
2**(1/12)
There are 12 semitones per octave,
\[ (2^{(1/12)})^{12} = 2\]
The major notes (white piano keys) are C, D, E, F, G, A, B, C5 (C5 is note C in the next octave.)
The MIDI note number for A above middle C is \(n=69\).
C, D, E, F, G, A, B, C5 = 60, 62, 64, 65, 67, 69, 71, 72
scale = [C, D, E, F, G, A, B, C5 ]
The minor notes (black piano keys) are the missing values: 61, 63, 66, 68, 70. We are going to ignore those in this simple demo.
Image('piano_keyboard_picture.jpg')
Let's play some sounds. First, A above middle C:
Fs = 22050 #samplerate
dur = 3 #note duration in seconds
t = linspace(0,dur,dur*Fs)
f = 440
Audio(cos(2*pi*f*t), rate=Fs)
Now let's play a scale:
BaseNote = 440
sound = []
for note in scale:
f = BaseNote * 2**((note-69)/12)
sound.append(cos(2*pi*f*t))
sound = concatenate(sound)
Audio(sound, rate=Fs)
Envelope Functions¶
The notes sound better if we apply an envelope function to soften the rise and fall. Professional envelopes use an ADSR (Attack, Delay, Sustain, Release) shape, but we will start simply.
The sine function starts at 0, rises to 1 (at pi/2), and returns to 0 (at pi).
def sine_env(t):
return sin(pi*t/t[-1]) #t[-1] = last value
plot(t,sine_env(t))
That works, but even better is to make it flatter. Here's a trick:
def flute_env(t):
return sin(pi*t/t[-1])**0.4
plot(t,flute_env(t),'b')
plot(t,sine_env(t), 'r')
All the code snippets have the same basic structure. sound=[] is an empty list. For each note in the song, the frequency is computed, a sinusoid of that frequency is created, and the sinusoid is appended to the list. When the for loop exits, the list of sounds is concatenated into a single longer sound. That sound is played with the Audio command.
BaseNote = 440
sound = []
for note in scale:
f = BaseNote * 2**((note-69)/12)
sinusoid = cos(2*pi*f*t)
sound.append(flute_env(t) * sinusoid)
sound = concatenate(sound)
Audio(sound, rate=Fs)
The sound from a stringed instrument (e.g., guitar or piano) rises quickly, then decreases slowly.
def guitar_env(t):
return (1-exp(-80*t))*exp(-8*t)
plot(t,guitar_env(t))
Here's an illustration of the guitar envelope times the sinusiod. The rise is very quick, almost impossible to see, but the decay is clear.
s = linspace(0,1,300)
plot(s,guitar_env(s)*cos(2*pi*30*s))
BaseNote = 440
dur = 0.5
sound = []
t = linspace(0,dur,dur*Fs)
for note in scale:
f = BaseNote * 2**((note-69)/12)
sinusoid = cos(2*pi*f*t)
sound.append(guitar_env(t) * sinusoid)
sound = concatenate(sound)
Audio(sound, rate=Fs)
Mary Had a Little Lamb¶
Let's play a song. Since the notes have different durations, we need to modify our code a bit. Each note has a frequency and duration. The song is a list of notes. When we create the song, we will loop over the notes.
mary = [ (E, 1/2), (D, 1/2), (C, 1/2), (D, 1/2), (E, 1/2), (E, 1/2),
(E, 1), (D, 1/2), (D, 1/2), (D, 1), (E, 1/2), (G, 1/2),
(G, 1), (E, 1/2), (D, 1/2), (C, 1/2), (D, 1/2), (E, 1/2),
(E, 1/2), (E, 1/2), (E, 1/2), (D, 1/2), (D, 1/2), (E, 1/2),
(D, 1/2), (C, 1)]
BaseNote = 440
sound = []
for note in mary:
fnum, dur = note
t = linspace(0,dur,dur*Fs)
f = BaseNote * 2**((fnum-69)/12)
sinusoid = cos(2*pi*f*t)
sound.append(guitar_env(t) * sinusoid)
gsound = concatenate(sound)
Audio(gsound, rate=Fs)
BaseNote = 880
sound = []
for note in mary:
fnum, dur = note
t = linspace(0,dur,dur*Fs)
f = BaseNote * 2**((fnum-69)/12)
sinusoid = cos(2*pi*f*t)
sound.append(flute_env(t) * sinusoid)
fsound = concatenate(sound)
Audio(fsound, rate=Fs)
#play guitar/piano and flute together
sound = gsound+0.3*fsound
Audio(sound, rate=Fs)
FM Modulation¶
Frequency modulation makes for an interesting sound.
\[X(t) = \cos\big(2\pi f t+k \sin(2\pi f_m t)\big)\]
Yamaha synthesizers were famous for their FM sound.
The snippet and plot below illustrate FM modulation. See how the frequency increases and decreases.
s = linspace(0,1,700)
plot(s, cos(2*pi*30*s + 5*sin(2*pi*5*s)) )
BaseNote = 220
fmsound = []
for note in mary:
fnum, dur = note
t = linspace(0,dur,dur*Fs)
f = BaseNote * 2**((fnum-69)/12)
FMsinusoid = cos(2*pi*f*t + 10*sin(4.1*pi*f*t))
fmsound.append( guitar_env(t) * FMsinusoid )
fmsound = concatenate(fmsound)
Audio(fmsound, rate=Fs)
sound = fmsound+0.5*fsound
Audio(sound, rate=Fs)
AM Modulation or Tremolo¶
Tremolo is an amplitude modulation of the note. The envelope has ripples. We must keep the tremolo frequency below the lower limit of hearing, 20 Hz.
\[ X(t) = \big(1+d\sin(2\pi f_r t)\big)\cos(2\pi f t)\]
\(d\) is the depth of the tremolo.
s = linspace(0,1,300)
tremolo = 1+0.6*sin(2*pi*5*s)
plot(s,tremolo*cos(2*pi*30*s))
BaseNote = 440
depth = 0.6
tremfreq = 8
sound = []
for note in mary:
fnum, dur = note
t = linspace(0,dur,dur*Fs)
tremolo = 1 + depth*sin(2*pi*tremfreq*t)
f = BaseNote * 2**((fnum-69)/12)
sinusoid = cos(2*pi*f*t)
sound.append(flute_env(t) * tremolo * sinusoid)
sound = concatenate(sound)
Audio(sound, rate=Fs)
Conclusion¶
Sinusoids are used for all sorts of things in Electrical & Computer Engineering.
We use complex algebra rather than trigonometry.
Sums of sinusoids at the same frequency result in a sinusoid at the same frequency.
Beat frequencies, modulation, frequency shifting result from combinations of sinusoids.
Other ideas: harmonics and Fourier series, phased arrays, noise cancellation.
Python is really handy. Start by using Python as a "super-calculator" and move on to programming concepts.