import numpy as np
import matplotlib.pyplot as plt

# 產(chǎn)生size點(diǎn)取樣的三角波，其周期為1
def triangle_wave(size):
  x = np.arange(0, 1, 1.0/size)
  y = np.where(x<0.5, x, 0)
  y = np.where(x>=0.5, 1-x, y)
  return x, y
 
def square_wave(size):
  x = np.arange(0, 1, 1.0/size)
  y = np.where(x<0.5, 1.0, 0)
  return x, y

fft_size = 256
 
# 計(jì)算三角波和其FFT
x, y = triangle_wave(fft_size)
fy = np.fft.fft(y) / fft_size

# 繪制三角波的FFT的前20項(xiàng)的振幅，由于不含下標(biāo)為偶數(shù)的值均為0， 因此取
# log之后無(wú)窮小，無(wú)法繪圖，用np.clip函數(shù)設(shè)置數(shù)組值的上下限，保證繪圖正確
plt.figure()
plt.plot(np.clip(20*np.log10(np.abs(fy[:20])), -120, 120), "o")
plt.xlabel("frequency bin")
plt.ylabel("power(dB)")
plt.title("FFT result of triangle wave")

# 取FFT計(jì)算的結(jié)果freqs中的前n項(xiàng)進(jìn)行合成，返回合成結(jié)果，計(jì)算loops個(gè)周期的波形
def fft_combine(freqs, n, loops=1):
  length = len(freqs) * loops
  data = np.zeros(length)
  index = loops * np.arange(0, length, 1.0) / length * (2 * np.pi)
  for k, p in enumerate(freqs[:n]):
    if k != 0: p *= 2 # 除去直流成分之外，其余的系數(shù)都*2
    data += np.real(p) * np.cos(k*index) # 余弦成分的系數(shù)為實(shí)數(shù)部
    data -= np.imag(p) * np.sin(k*index) # 正弦成分的系數(shù)為負(fù)的虛數(shù)部
  return index, data

# 繪制原始的三角波和用正弦波逐級(jí)合成的結(jié)果，使用取樣點(diǎn)為x軸坐標(biāo)
plt.figure()
plt.plot(y, label="original triangle", linewidth=2)
for i in [0,1,3,5,7,9]:
  index, data = fft_combine(fy, i+1, 2) # 計(jì)算兩個(gè)周期的合成波形
  plt.plot(data, label = "N=%s" % i)
plt.legend()
plt.title("partial Fourier series of triangle wave")
plt.show()