%pylab inline
from __future__ import print_function
from __future__ import division
Populating the interactive namespace from numpy and matplotlib
x = linspace(0, 2*pi, 500)
plot(x, sin(x))
[<matplotlib.lines.Line2D at 0x7f87a7967890>]
http://en.wikipedia.org/wiki/Trigonometric_functions
The Sine function relates the length of the opposite side to the angle of a right angle triangle:
$$\sin(\theta)= \frac{opposite}{hypotenuse}$$The Cosine function is the sine function with a 90 degree ($\frac{\pi}{2}$) phase shift. It is related to the adjacent side in a right angle triangle:
$$\cos(\theta)= \frac{adjacent}{hypotenuse}$$A phase vector with constant amplitude and angular velocity.
An angle (for a constant line of length $A$) changes at a constant speed.
http://en.wikipedia.org/wiki/Phasor
$$A\cdot \cos(\omega t + \theta)$$Now discrete:
theta = linspace(0,2*pi, 10, endpoint=False)
print(theta)
[ 0. 0.62831853 1.25663706 1.88495559 2.51327412 3.14159265 3.76991118 4.39822972 5.02654825 5.65486678]
theta = linspace(0,2*pi, 10, endpoint=False)
r = ones_like(theta) # ones(theta.shape)
polar(theta, r, 'o')
for ang in theta:
annotate('{:.2f}$\pi$'.format(ang/pi), xy=(ang,1))
#xticks(())
#yticks(());
#grid()
r
array([ 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.])
The separation in angles is constant, so if the separation between "samples" a constant time, then you can say the angular velocity is constant.
phs1 = linspace(0, 5 * 2*pi, 300)
y1 = sin(phs1)
y2 = 0.7 * sin(phs1 + 0.7*pi)
plot(y1)
plot(y2);
plot(y1 + y2);
plot(y1)
plot(y2);
plot(y1 + y2);
phs1 = linspace(0, 5 * 2*pi, 300)
y1 = sin(phs1)
y2 = 0.7 * sin(phs1 + 1.2*pi)
plot(y1)
plot(y2);
plot(y1 + y2);
phs1 = linspace(0, 5 * 2*pi, 300)
y1 = sin(phs1)
y2 = 0.7 * sin(phs1 + 0.5*pi)
plot(y1)
plot(y2);
plot(y1 + y2);
phs1 = linspace(0, 5 * 2*pi, 300)
y1 = sin(phs1)
y2 = 0.7 * sin(phs1 + 1*pi)
plot(y1)
plot(y2);
plot(y1 + y2);
Adding two sinusoids of the same frequency will always result in the same frequency, no matter the phase.
The resulting amplitude of the sum depends on the original sinusoids' amplitude and their phase relationships
You can easily calculate this with phasors. Adding sinusoids is the same as adding the rotating phasor vectors.
You can see:
But, the most interesting implication is that by adding a single phasor with initial phase of 0 with another with a 90 degree shift, a phasor with any phase and amplitude can be described. i.e.:
$$ A \cdot \sin(\omega t + \phi) = [A \cdot \cos(\phi) ]\sin(\omega t) + [A \cdot \sin(\phi)]\cos(\omega t)$$This comes from the trigonometric identity:
$$\sin(a + b) = \sin(a)\cos(b) + \cos(a)\sin(b)$$What if we consider the x-y plane, as the complex plane, where the x-axis describes the real part of a number and the y-axis the imaginary part? Then we can describe any phasor as:
$$ a \cdot \cos(\omega t) + b \cdot j \cdot \sin(\omega t)$$http://en.wikipedia.org/wiki/Euler%27s_formula
$$\ \ \ \ \ \ \ e^{jx} = \cos(x) + j\sin(x)$$So Euler's formula describes a phasor in two different ways!
And:
$$\cos(\omega t) = \mathrm{Re}[ e^{j\omega t}] $$$$\sin(\omega t) = \mathrm{Im}[ e^{j\omega t}] $$$$ \ \ \ \ e^{j\pi} + 1 = 0 $$amp_harm_rel = lambda x: 1.0/x
print(amp_harm_rel)
amp_harm_rel(5)
<function <lambda> at 0x7f87a7a742a8>
0.2
Equivalent to:
def amp_harp_rel_func(x):
return 1.0/x
print(amp_harp_rel_func)
amp_harp_rel_func(5)
<function amp_harp_rel_func at 0x7f87a5ed1848>
0.2
amp_harm_rel = lambda x: 1.0/x
harmonics = arange(30) + 1
print(harmonics)
[ 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30]
amps = amp_harm_rel(harmonics)
plot(amps, 'o')
[<matplotlib.lines.Line2D at 0x7f87a5d36650>]
zip(harmonics,amps)
[(1, 1.0), (2, 0.5), (3, 0.33333333333333331), (4, 0.25), (5, 0.20000000000000001), (6, 0.16666666666666666), (7, 0.14285714285714285), (8, 0.125), (9, 0.1111111111111111), (10, 0.10000000000000001), (11, 0.090909090909090912), (12, 0.083333333333333329), (13, 0.076923076923076927), (14, 0.071428571428571425), (15, 0.066666666666666666), (16, 0.0625), (17, 0.058823529411764705), (18, 0.055555555555555552), (19, 0.052631578947368418), (20, 0.050000000000000003), (21, 0.047619047619047616), (22, 0.045454545454545456), (23, 0.043478260869565216), (24, 0.041666666666666664), (25, 0.040000000000000001), (26, 0.038461538461538464), (27, 0.037037037037037035), (28, 0.035714285714285712), (29, 0.034482758620689655), (30, 0.033333333333333333)]
zip(harmonics,arange(20))
[(1, 0), (2, 1), (3, 2), (4, 3), (5, 4), (6, 5), (7, 6), (8, 7), (9, 8), (10, 9), (11, 10), (12, 11), (13, 12), (14, 13), (15, 14), (16, 15), (17, 16), (18, 17), (19, 18), (20, 19)]
x = linspace(0, 2*pi, 500)
out = zeros_like(x)
for harm, amp in zip(harmonics,amps):
out += amp * sin(x*harm)
plot(out)
[<matplotlib.lines.Line2D at 0x7f87a5df77d0>]
def add_harmonics(num, amp_harm_rel, numpoints=500):
x = linspace(0, 2*pi, numpoints)
out = zeros_like(x)
harmonics = arange(num) + 1
amps = amp_harm_rel(harmonics)
for harm, amp in zip(harmonics,amps):
out += amp * sin(x*harm)
return out
amp_harm_rel = lambda x: 1.0/x
for i in range(20):
plot(add_harmonics(i + 1, amp_harm_rel))
amp_harm_rel = lambda x: 1.0/x
for i in range(20):
plot(add_harmonics(i + 10, amp_harm_rel))
plot(add_harmonics(300, amp_harm_rel))
[<matplotlib.lines.Line2D at 0x7f87a59f02d0>]
Adding sine wave and its harmonics at amplitude of $\frac{1}{N}$ produces a saw wave.
We also saw that adding sine waves with the same amplitude results in a bi-polar impulse:
def amp_harm_rel(x):
return ones_like(x)
for i in range(20):
plot(add_harmonics(i + 10, amp_harm_rel))
Adding only the odd harmonics witha amplitude of $\frac{1}{N}$ produces a square wave:
amp_harm_rel = (lambda x: ceil(x/2.0 - (x/2).astype(int)) * 1.0/x)
for i in range(20):
plot(add_harmonics(i + 1, amp_harm_rel))
amp_harm_rel = (lambda x: ceil(x/2.0 - (x/2).astype(int)) * 1.0/x)
for i in range(20):
plot(add_harmonics(i + 20, amp_harm_rel))
These are all examples of Fourier Series, which allow calculating a continuous periodic function through the sum of simple sinusoids.
A Fourier series is described as:
$$s_N(x) = \frac{a_0}{2} + \sum_{n=1}^N A_n\cdot \sin(\tfrac{2\pi nx}{P}+\phi_n)$$Fourier's theorem
Any periodic function can be written as a Fourier series. i.e. any periodic function can be written as a sum of harmonic phasors.
$$S(f) = \int_{-\infty}^{\infty} s(t) \cdot e^{- i 2\pi f t} dt$$Which is multiplying the input function by a phasor of the "input" frequency.
When working with discrete data, you need to use the Discrete Fourier Transform (DFT):
$$X_k = \sum_{n=0}^{N-1} x_n \cdot e^{-i 2 \pi k n / N}$$phs = linspace(0, 10* 2 * pi, 512, endpoint=False)
x = cos(phs)
plot(phs, x)
[<matplotlib.lines.Line2D at 0x7f87a5732050>]
phasor_phs = linspace(0, 2 * pi, 512, endpoint=False)
plot(x)
plot(sin(phasor_phs))
plot(cos(phasor_phs))
[<matplotlib.lines.Line2D at 0x7f87a5a85d90>]
60*44100
2646000
plot(x*sin(phasor_phs))
plot(x*cos(phasor_phs))
sum(x*sin(phasor_phs)), sum(x*cos(phasor_phs))
(-2.7474550412520671e-14, -2.2537527399890678e-14)
plot(x*sin(2*phasor_phs))
plot(x*cos(2*phasor_phs))
sum(x*sin(2*phasor_phs)), sum(x*cos(2*phasor_phs))
(-3.0253577421035516e-15, -1.8873791418627661e-14)
plot(x*sin(5*phasor_phs))
plot(x*cos(5*phasor_phs))
sum(x*sin(5*phasor_phs)), sum(x*cos(5*phasor_phs))
(-1.865174681370263e-14, -1.7652546091539989e-14)
dft = []
for k in range(len(x)):
fft_bin = complex(sum(x*cos(k * phasor_phs)), sum(-x*sin(k * phasor_phs)))
dft.append(fft_bin)
subplot(121)
plot(real(dft))
title('Real part')
subplot(122)
plot(imag(dft))
title('Imaginary part')
gcf().set_figwidth(10)
plot(real(dft))
xlim((0, 20))
(0, 20)
plot(real(dft))
xlim((500,512))
(500, 512)
phs_alias = linspace(0, 35* 2*pi, 40, endpoint=False)
plot(sin(phs_alias))
[<matplotlib.lines.Line2D at 0x7f87a562a310>]
The hypotenuse of the complex number in the complex plane is the magnitude of the component and the angle is the phase. If you plot all hypotenuses, you get what's called the Magnitude spectrum:
plot(hypot(real(dft), imag(dft)))
[<matplotlib.lines.Line2D at 0x7f87a553ced0>]
plot(hypot(real(dft), imag(dft)))
xlim((0, 50))
plot(real(fft.fft(x)))
[<matplotlib.lines.Line2D at 0x7f87a5cacbd0>]
plot(imag(fft.fft(x)))
[<matplotlib.lines.Line2D at 0x7f87a5bae950>]
plot(abs(fft.fft(x)))
title('The magnitude spectrum')
<matplotlib.text.Text at 0x7f87a58ae090>
plot(angle(fft.fft(x)))
title('The phase spectrum')
<matplotlib.text.Text at 0x7f87a5d28f10>
plot(angle(fft.fft(x)))
xlim((0, 20))
(0, 20)
plot(angle(fft.fft(sin(phs))))
xlim((0, 20))
(0, 20)
By: Andrés Cabrera mantaraya36@gmail.com For MAT course MAT 201A at UCSB
This ipython notebook is licensed under the CC-BY-NC-SA license: http://creativecommons.org/licenses/by-nc-sa/4.0/