%matplotlib inline
import matplotlib.pyplot as plt
import matplotlib.image as mpimg
import numpy as np
import scipy.io.wavfile as wavfile
from scipy import integrate
plt.rcParams['figure.figsize'] = (15, 10)
# sampling_rate, signal = wavfile.read("flute.wav")
# x = np.arange(len(signal))
# testing with Andreas' toy signal
y = lambda x: (0.5 if x < 0.25 or x > 0.75 else -0.5 )+ 0.5
x = np.arange(0,1.0,0.001)
signal = np.array([y(i) for i in x])
x.shape, signal.shape
((1000,), (1000,))
plt.plot(signal)
[<matplotlib.lines.Line2D at 0x123f75510>]
# Our lowest frequency should be the length of our signal.
# We assume that the signal is periodic with the first period being the length of our signal.
T = len(signal)
f0 = 1.0/T
def a_coeff(x,y,i,f0,c):
return f0*2.0*integrate.trapz((y-c)*np.sin(2*np.pi*x*f0*i),x)
def b_coeff(x,y,i,f0,c):
return f0*2.0*integrate.trapz((y-c)*np.cos(2*np.pi*x*f0*i),x)
def c_coeff(x,y):
return integrate.trapz(y,x)
c = c_coeff(x,signal)
a_s = []
b_s = []
for i in xrange(1,100):
a_s.append(a_coeff(x,signal,i,f0,c))
b_s.append(b_coeff(x,signal,i,f0,c))
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 xrange(0,i_max):
s += a_s[i]*np.sin(2*np.pi*x*f0*(i+1))+b_s[i]*np.cos(2*np.pi*x*f0*(i+1))
return s
approximation = [fourier_value(i,a_s,b_s,c,f0,50) for i in x]
plt.plot(approximation)
[<matplotlib.lines.Line2D at 0x12409c6d0>]