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
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
#!/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()