%pylab inline
from matplotlib.pylab import *
Populating the interactive namespace from numpy and matplotlib
from numpy import *
import math
x = arange(0,2*math.pi,0.01)
plot(x/math.pi/2.0,sin(x));
plot(x/math.pi/2.0,sin(x*2));
plot(x/math.pi/2.0,sin(3*x)*0.5+sin(x*2)+sin(x)*2);
axvline(0.5,ls = '--')
xlabel("t / $f_0$")
<matplotlib.text.Text at 0x105e2b4d0>
from numpy import *
import math
x = arange(0,2*math.pi,0.01)
plot(x/math.pi/2.0,cos(x));
plot(x/math.pi/2.0,cos(x*2));
plot(x/math.pi/2.0,cos(3*x)*0.5+cos(x*2)+cos(x)*2);
axvline(0.5,ls = '--')
xlabel("t / $f_0$")
<matplotlib.text.Text at 0x105e60c50>
from numpy import *
import math
x = arange(0,2*math.pi,0.001)
plot(x/math.pi/2.0,cos(x)*cos(x));
axvline(0.5,ls = '--')
xlabel("t / $f_0$")
sum(cos(x)*cos(x))
3142.4073464260332
sy = lambda x: (0.5 if x < 0.25 or x > 0.75 else -0.5 )+ 0.5
sx = arange(0,1.0,0.001)
plot(sx,[sy(x) for x in sx])
#ylim(-0.6,0.6)
signal = sx,[sy(x) for x in sx]
f0 = 1.0/signal[0][-1]
from scipy.integrate import trapz
def a_coeff(signal,i,f0,c):
x,y = signal
return f0*2.0*trapz((y-c)*sin(2*math.pi*x*f0*i),x)
def b_coeff(signal,i,f0,c):
x,y = signal
return f0*2.0*trapz((y-c)*cos(2*math.pi*x*f0*i),x)
def c_coeff(signal):
x,y = signal
return trapz(y,x)
c = c_coeff(signal)
a_s = []
b_s = []
for i in range(1,100):
a_s.append(a_coeff(signal,i,f0,c))
b_s.append(b_coeff(signal,i,f0,c))
c,a_s[:10],b_s[:10]
(0.498, [0.0020019797286802026, 1.8887031483585555e-05, -0.0020018015450786646, -3.7772381950862726e-05, 0.0020014451937346088, 5.6654370536769179e-05, -0.0020009107063646159, -7.5531316674086157e-05, 0.0020001981305400539, 9.4401540244977599e-05], [0.63661059107062334, 0.0030029188598676843, -0.2121790477406521, -0.0030026664386361398, 0.12727805220031627, 0.0030022457638306924, -0.090881424480230483, -0.0030016568763193601, 0.070652923858416386, 0.0030008998333135536])
def fourier_value(x,a_s,b_s,c,f0,i_max = None ):
s = c
i_max = i_max or len(a_s)
for i in range(0,i_max):
s += a_s[i]*sin(2*math.pi*x*f0*(i+1))+b_s[i]*cos(2*math.pi*x*f0*(i+1))
return s
fourier_value(0.0,a_s,b_s,c,f0,None)
0.99816089257685292
plot(signal[0],[fourier_value(x,a_s,b_s,c,f0,50) for x in signal[0]])
plot(signal[0],signal[1])
#ylim(-0.1,1.1)
[<matplotlib.lines.Line2D at 0x10604da10>]