In [ ]:
In [33]:
# Compute the Discrete Fourier Transform of the sequence of numbers l
def dft(l):
    N=len(l)
    output=[]
    k=0
    while k<=N-1:
        n=0
        sum=0
        while n<=N-1:
           sum=sum+l[n]*e^(-I*2*pi/N*k*n)
           n=n+1
        output.append(sum)
        k=k+1
    return output
In [2]:
# Compute the DFT of the sequence a, but just print the norms of 
# the complex numbers
def norm_dft(a):
    u=dft(a)
    output=[]
    for v in u:
        x=N(v.abs())
        if x<0.001:
            x=0
        output.append(x)
    return output
In [56]:
f=piecewise([((0,2), 0), ((2,3), 1), ((3,5),0), ((5,6),1), ((6,8),0), ((8,9),1),((9,11),0)])
In [57]:
plot(f,(0,11))  # This is a function defined on [0,11], but with period 3
Out[57]:
In [50]:
plot(f,(0,11))+plot(fourier_approx(f,15,11,x),(x,0,11),color='red')+plot(fourier_approx(f,4,11,x),(x,0,11),color='green')
Out[50]:
In [3]:
# This function gives the coefficients in the Fourier series for f.
# Prints out [k:a_k,b_k,f_k] where a_k is the coefficient of cosine,
# b_k is the coefficient of sine, and f_k is the frequency
#
def fourier_exp(f,n,length):    #compute series from k=0 to k=n
    i=1
    L=length/2
    a0=f.fourier_series_cosine_coefficient(0)/2
    output=[[0,':',a0,0,0]]
    while i<=n:
        a=N(f.fourier_series_cosine_coefficient(i))
        if abs(a)<0.0001:
            a=0
        b=N(f.fourier_series_sine_coefficient(i))
        if abs(b)<0.0001:
            b=0
        output.append([i,':',a,b,i/L*pi])
        i=i+1
    return output
In [47]:
fourier_exp(f,15,11)
Out[47]:
[[0, ':', 3/11, 0, 0],
 [1, ':', -0.128306329010999, 0, 2/11*pi],
 [2, ':', -0.158149502028013, 0, 4/11*pi],
 [3, ':', -0.293619439869413, 0, 6/11*pi],
 [4, ':', 0.388353062244215, 0, 8/11*pi],
 [5, ':', 0.0390335710237439, 0, 10/11*pi],
 [6, ':', -0.0325279758531199, 0, 12/11*pi],
 [7, ':', -0.221916035568123, 0, 14/11*pi],
 [8, ':', 0.110107289951030, 0, 16/11*pi],
 [9, ':', 0.0351443337840029, 0, 18/11*pi],
 [10, ':', 0.0128306329010999, 0, 20/11*pi],
 [11, ':', 0, 0, 2*pi],
 [12, ':', -0.0106921940842499, 0, 24/11*pi],
 [13, ':', -0.0243306926196943, 0, 26/11*pi],
 [14, ':', -0.0629184514005885, 0, 28/11*pi],
 [15, ':', 0.103560816598457, 0, 30/11*pi]]
In [ ]:
In [84]:
# The gives the Fourier approximation function defined using a sum
# from k=0 to k=n
# Here "length" is the length of the original interval on which f
# is defined.
#
def fourier_approx(f,n,length,x):
    list=fourier_exp(f,n,length)
    temp=list[0][2]
    i=1
    while i<=n:
        freq=list[i][4]
        A=list[i][2]
        B=list[i][3]
        temp=temp+A*cos(freq*x)+B*sin(freq*x)
        i=i+1
    return temp
In [59]:
f=piecewise([ ((0,1),0),((1,2),1 ) ])
In [61]:
plot(f,(0,2))
Out[61]:
In [68]:
plot(f,(0,2))+plot(fourier_approx(f,1,2,x),(x,0,2),color='red')
Out[68]:
In [65]:
plot(f,(0,2))+plot(fourier_approx(f,5,2,x),(x,0,2),color='red')
Out[65]:
In [69]:
f=[0,0,1,0,0,1,0,0,1,0,0]   #Here N=11 and s=3 
norm_dft(f)                 #Get spikes in F near 11/3 and 2*(11/3)
Out[69]:
[3.00000000000000,
 0.715370323453430,
 0.918985947228995,
 1.83083002600377,
 2.68250706566236,
 0.309721467890568,
 0.309721467890572,
 2.68250706566236,
 1.83083002600377,
 0.918985947228996,
 0.715370323453425]
In [72]:
f=[0,0,1,0,0,1,0,0,1]
f=f+f+f                  # Here N=56 and s=3.  Get spikes in F near
f=f+f+[0,0]              # 56/3 and 2*(56/3)
norm_dft(f)
Out[72]:
[18.0000000000000,
 0.668419798944521,
 0.673736120718789,
 0.682791329343547,
 0.695895486700944,
 0.713521983562587,
 0.736355241089115,
 0.765366864730180,
 0.801937735804843,
 0.848058351891489,
 0.906669817410423,
 0.982273197446115,
 1.08208834612854,
 1.21844027553862,
 1.41421356237310,
 1.71723086396957,
 2.24697960371747,
 3.40954570983426,
 8.04691719426355,
 15.1010780357509,
 3.51351878929970,
 1.84775906502257,
 1.17190265282955,
 0.798431468126536,
 0.554958132087373,
 0.377365428858003,
 0.235750464926204,
 0.113569091696766,
 0,
 0.113569091696802,
 0.235750464926220,
 0.377365428857995,
 0.554958132087371,
 0.798431468126585,
 1.17190265282957,
 1.84775906502257,
 3.51351878929973,
 15.1010780357508,
 8.04691719426358,
 3.40954570983431,
 2.24697960371745,
 1.71723086396955,
 1.41421356237310,
 1.21844027553866,
 1.08208834612850,
 0.982273197446101,
 0.906669817410410,
 0.848058351891429,
 0.801937735804848,
 0.765366864730180,
 0.736355241089097,
 0.713521983562637,
 0.695895486700941,
 0.682791329343491,
 0.673736120718728,
 0.668419798944479]
In [73]:
list_plot(f)+list_plot(norm_dft(f),color='red')
Out[73]:
In [90]:
f=piecewise([((0,2), 0), ((2,3), 1), ((3,5),0), ((5,6),1), ((6,8),0), ((8,9),1),((9,11),0)])
fourier_exp(f,25,11)    #N=11, s=3: note the biggest contributions are
                       #around k=11/3, 2*(11/3), 3*(11/3), 4*(11/3) and so on
Out[90]:
[[0, ':', 3/11, 0, 0],
 [1, ':', -0.128306329010999, 0, 2/11*pi],
 [2, ':', -0.158149502028013, 0, 4/11*pi],
 [3, ':', -0.293619439869413, 0, 6/11*pi],
 [4, ':', 0.388353062244215, 0, 8/11*pi],
 [5, ':', 0.0390335710237439, 0, 10/11*pi],
 [6, ':', -0.0325279758531199, 0, 12/11*pi],
 [7, ':', -0.221916035568123, 0, 14/11*pi],
 [8, ':', 0.110107289951030, 0, 16/11*pi],
 [9, ':', 0.0351443337840029, 0, 18/11*pi],
 [10, ':', 0.0128306329010999, 0, 20/11*pi],
 [11, ':', 0, 0, 2*pi],
 [12, ':', -0.0106921940842499, 0, 24/11*pi],
 [13, ':', -0.0243306926196943, 0, 26/11*pi],
 [14, ':', -0.0629184514005885, 0, 28/11*pi],
 [15, ':', 0.103560816598457, 0, 30/11*pi],
 [16, ':', 0.0121979909449200, 0, 32/11*pi],
 [17, ':', -0.0114804620658070, 0, 34/11*pi],
 [18, ':', -0.0863006804987144, 0, 36/11*pi],
 [19, ':', 0.0463609641899073, 0, 38/11*pi],
 [20, ':', 0.0158149502028013, 0, 40/11*pi],
 [21, ':', 0.00610982519099997, 0, 42/11*pi],
 [22, ':', 0, 0, 4*pi],
 [23, ':', -0.00557853604395649, 0, 46/11*pi],
 [24, ':', -0.0131791251690011, 0, 48/11*pi],
 [25, ':', -0.0352343327843296, 0, 50/11*pi]]
In [93]:
# Here we show f in blue, the Fourier approx in red, and the 
#"largest coefficient" terms of the Fourier approximation in green
plot(f,(0,11))+plot(fourier_approx(f,25,11,x),(0,11),color='red')+plot(3/11-0.29*cos(6/11*pi*x)+0.388*cos(8/11*pi*x)-0.22*cos(14/11*pi*x)+0.11*cos(16/11*pi*x),(0,11),color='green')
Out[93]:
In [37]:
# Example of a sum of cosine and sine functions with different periods
plot(cos(3*x)+sin(5*x)+6*cos(x),(0,20))
Out[37]:
In [ ]:
In [ ]: