#!/usr/bin/env python
import scipy.integrate as si
import pylab as s
import threading as t
def signal_func ( x ):
return 2 * s.cos ( 2 * x ) + .5 * s.cos ( 1.5 * x ) - s.sin ( .5 * x )
def f_analysis ( x, func, order=None ):
coeff_0, coeff_1, omega = [], [], []
limit = order
if limit == None:
limit = len(x)
for n in range ( limit ):
omega.append ( 2 * s.pi * n / (x[-1]-x[0]) )
coeff_0.append ( 2./ (x[-1]-x[0]) * si.romb( func(x) \
* s.cos ( omega[n] * x ), abs(x[1]-x[0]) ) )
coeff_1.append ( 2./ (x[-1]-x[0]) * si.romb( func(x) \
* s.sin ( omega[n] * x ), abs(x[1]-x[0]) ) )
return omega, coeff_0, coeff_1
def f_row_print ( omega, coeff0, coeff1, order = 5):
print 'Coefficients of Fourier Row:'
for i in range ( max(len(omega),order) ):
if abs(coeff0[i]) > 0.009:
print 'a[%i] = %4.2f * cos ( %4.2f * x )' % (i, coeff0[i], omega[i])
if abs(coeff1[i]) > 0.009:
print 'b[%i] = %4.2f * sin ( %4.2f * x )' % (i, coeff1[i], omega[i])
def f_synthesis ( x, omega, coeff_0, coeff_1 ):
y = s.zeros(len(x), dtype='Float64')
for i in range ( len(omega) ):
y += coeff_0[i] * s.cos ( omega[i] * x ) \
+ coeff_1[i] * s.sin ( omega[i] * x )
return y
if __name__ == '__main__':
x = s.linspace ( -2*s.pi, 2*s.pi, 2**10 + 1 )
y = signal_func ( x )
s.figure()
s.subplot(311)
s.plot ( x, signal_func(x), 'r-' )
f_row = f_analysis ( x, signal_func, 100 )
f_row_print ( f_row[0], f_row[1], f_row[2] )
s.subplot(312)
s.plot ( f_row[0], f_row[1], 'go', f_row[0], f_row[2], 'bo' )
F_row = f_synthesis ( x, f_row[0], f_row[1], f_row[2] )
s.subplot(313)
s.plot ( x, F_row )
s.show()