我們經常使用傅立葉變換來計算數字訊號的頻譜，進而分析數字訊號，離散時間傅立葉變換的公式為：
$$

可是自己動手實現一遍才是最好的學習。

在數字分析裡面，傅立葉變換預設等時間間隔取樣，不需要時間序列，只需要訊號陣列即可分析。

分析過程如下：

• 對於含有 n 個樣本值的數字訊號序列，根據奈奎斯特取樣定律，包含的週期數最大為 n/2,週期數為 0 代表直流分量。所以，當週期數表示為離散的 0,1,2,3…n/2 ，總的數目為 n/2+1個
• 傅立葉變換之後的結果為複數， 下標為 k 的複數 a+b*j 表示時域訊號中週期為 N/k 個取樣值的正弦波和餘弦波的成分的多少， 其中 a 表示 cos 波形的成分， b 表示 sin 波形的成分
• 首先產生一個長度為 n,一倍週期的 $e^{-jwn} $ （即為 $cos(wn)-jsin(wn) $ ）波樣本序列.
• 將數字訊號序列中的每一個樣本與 1 倍週期的樣本波形序列相乘，得到 n 個乘積，將 n 個乘積相加，放入 f[1] 中。
• 再產生一個長度為 n,兩倍週期的 $e^{-jwn} $ （即為 $cos(wn)-jsin(wn) $ ）波樣本序列,再將數字訊號序列中的每一個樣本與 2 倍週期的樣本波形序列相乘，得到 n 個乘積，將 n 個乘積相加，放入 f[2] 中。依次重複。
• 對於 0 倍週期，即直流分量來說，可以認為產生的是 0 倍週期的樣本波形，重複操作，放入 f[0] 即可。
• 這樣就得到了數字訊號序列的傅立葉變換

使用方法：

從以上過程得到數字序列的傅立葉變換之後，如果想要得到真正頻譜，還需要做處理：

• 計算出的每一個頻率下的幅值需要除以時間序列的長度，類似求平均的過程
• 每一個頻率下的幅值是一個複數，需要對它求模，而且因為在負頻率處也有值，所以需要對於實訊號需要乘 2
• 頻率的序列為 0 到取樣率的一半，長度為 n/2+1
  

完整程式：

# 離散時間傅立葉變換的 python 實現
import numpy as np
import math
import pylab as pl
import scipy.signal as signal
import matplotlib.pyplot as plt

sampling_rate=1000
t1=np.arange(0,10.0,1.0/sampling_rate)
x1 =np.sin(15*np.pi*t1)

# 傅立葉變換
def fft1(xx):
#   t=np.arange(0,s)
  t=np.linspace(0,1.0,len(xx))
  f = np.arange(len(xx)/2+1,dtype=complex)
  for index in range(len(f)):
    f[index]=complex(np.sum(np.cos(2*np.pi*index*t)*xx),-np.sum(np.sin(2*np.pi*index*t)*xx))
  return f

# len(x1)

xf=fft1(x1)/len(x1)
freqs = np.linspace(0,sampling_rate/2,len(x1)/2+1)
plt.figure(figsize=(16,4))
plt.plot(freqs,2*np.abs(xf),'r--')

plt.xlabel("Frequency(Hz)")
plt.ylabel("Amplitude($m$)")
plt.title("Amplitude-Frequency curve")

plt.show()

plt.figure(figsize=(16,'r--')

plt.xlabel("Frequency(Hz)")
plt.ylabel("Amplitude($m$)")
plt.title("Amplitude-Frequency curve")
plt.xlim(0,20)
plt.show()

此處實現的是傳統的傅立葉變換，這種方法實際已經不用了，現在使用快速傅立葉變換，其實兩種是等價的，但是快速傅立葉變換時間複雜度要小很多。

以上就是本文的全部內容，希望對大家的學習有所幫助，也希望大家多多支援我們。