1. 程式人生 > 程式設計 >Python FFT合成波形的例項

Python FFT合成波形的例項

使用Python numpy模組帶的FFT函式合成矩形波和方波,增加對離散傅立葉變換的理解。

匯入模組

import numpy as np
import matplotlib.pyplot as plt


分別是產生一個週期的方波和三角波程式

# 產生size點取樣的三角波,其週期為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.0,0)
  return x,y


其中np.where函式第二個值是if,第三個是else

下面程式可以計算對應的頻譜,取樣點數取為2的n次冪是為了便於FFT計算

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


下面對計算的頻譜進行視覺化,頻率對應的強度使用工程上常用的分貝dp來表示

# 繪製三角波的FFT的前20項的振幅,由於不含下標為偶數的值均為0, 因此取
# log之後無窮小,無法繪圖,用np.clip函式設定陣列值的上下限,保證繪圖正確
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計算的結果freqs中的前n項進行合成,返回合成結果,計算loops個週期的波形
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 # 除去直流成分之外,其餘的係數都*2
    data += np.real(p) * np.cos(k*index) # 餘弦成分的係數為實數部
    data -= np.imag(p) * np.sin(k*index) # 正弦成分的係數為負的虛數部
  return index,data

其中index代表頻譜空間的取樣點

畫出合成訊號,x座標使用預設的整數表示即可

# 繪製原始的三角波和用正弦波逐級合成的結果,使用取樣點為x軸座標
plt.figure()
plt.plot(y,label="original triangle",linewidth=2)
for i in [0,3,5,7,9]:
  index,data = fft_combine(fy,i+1,2) # 計算兩個週期的合成波形
  plt.plot(data,label = "N=%s" % i)
plt.legend()
plt.title("partial Fourier series of triangle wave")
plt.show()

以上這篇Python FFT合成波形的例項就是小編分享給大家的全部內容了,希望能給大家一個參考,也希望大家多多支援我們。