今天再來做一些昨天方波傅立葉級數分析的延伸。

我們想要看到傅立葉級數分析下，正弦/餘弦波是怎麼樣組合成方波。

於是改寫程式碼如下：

import numpy as np
import matplotlib.pyplot as plt

def square_wave(T, t, D): #T週期, t時間點, D方波在高波的時間比例
    tinT = np.mod(t, T) #找出t在週期中是哪個階段，mod為餘數函數，取t除以Ｔ之後的餘數
    if tinT <= D*T:
        return 1 #若該時間在週期中處於高波時間比例的範圍內，則得出1，1為高波震幅
    else:
        return -1 #不在高波時間比例的範圍內則得出-1，-1為低頻的震幅。

#製圖精度
N = 1000

#定義週期
T = 10

#時間的最大值
t_max = T

#定義方波的高波比例
D = 0.5

#時間範圍
t_list = [t_max*i/N for i in range(N)]
#空間範圍
field_list = [square_wave(T, t_list[i], D) for i in range(N)]

#兩個函數做內積（兩個函數相乘對空間做積分
#這裡直接將兩個函數的內容定義為陣列array去做對位相乘
def inner_product(list1, list2, dx):
    #先將兩筆函數內容list定義為array型式
    array1 = np.array(list1)
    array2 = np.array(list2)
    product = array1*array2*dx
    return product.sum() #將相乘結果加總

#設定cos傅立葉級數
def Fourier_series_cos(field_list, n):
    cos_list = [np.cos(n*f0*t_list[i]) for i in range(N)]
    return inner_product(cos_list, field_list, T/N)

#設定sin傅立葉級數
def Fourier_series_sin(field_list, n):
    sin_list = [np.sin(n*f0*t_list[i]) for i in range(N)]
    return inner_product(sin_list, field_list, T/N)

#設定基礎頻率
f0 = 2*np.pi/T

#定義傅立葉級數序號，可自行帶入數值，這裏示範為10
n_max = 10
n_list = [i for i in range(n_max)]

#定義傅立葉級數轉換結果cos list
Fourier_series_list_cos = [Fourier_series_cos(field_list, n_list[i]) for i in range(n_max)]

#定義傅立葉級數轉換結果sin list
Fourier_series_list_sin = [Fourier_series_sin(field_list, n_list[i]) for i in range(n_max)]


#建立繪圖迴圈
plt.subplot(1,2,1)
plt.plot(t_list, field_list)
sum_array = np.zeros(N) #定義陣列及為陣列
for n in range(n_max):
    sin_list = [Fourier_series_list_sin[n]*np.sin(n*f0*t_list[i]) for i in range(N)]
    sin_array = np.array(sin_list) #將sin_list做成array
    sum_array += sin_array #加總arrays
    plt.plot(t_list, sin_list)

plt.plot(t_list, sum_array, color = 'deepskyblue')

plt.subplot(1,2,2)
plt.bar(n_list, Fourier_series_list_sin)
plt.show()

跑出的圖如下：

程式碼變更內容如下：

1. 新增n_max，這樣只需更改一處的傅立葉級數序號即可
2. 方波off狀態震幅改為-1
3. 建立繪圖迴圈，以畫出級數範圍內n的組合圖，這裏以sin為例
4. 改成條狀圖視覺效果更佳

上圖中可以看到，左圖天空藍代表sin arrays加總的結果，而這條天空藍的波由其他sin的波組合而成。

另外，如果須得出傅立葉級數轉換內積和數值，可以加入以下程式碼：

#定義arrays
cos1_list = [np.cos(f0*t_list[i]) for i in range(N)]
cos2_list = [np.cos(2*f0*t_list[i]) for i in range(N)]#兩倍波頻 n =2
sin1_list = [np.sin(f0*t_list[i]) for i in range(N)]
sin2_list = [np.sin(2*f0*t_list[i]) for i in range(N)]#兩倍波頻 n =2

#assign 內積函數結果
value = inner_product(cos1_list, cos2_list, T/N)

print(value)

好的，那麼傅立葉級數轉換，以方波為例就討論到這。
明天會開始新的篇章唷！敬請期待！