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

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

在數字分析里面，傅里葉變換默認等時間間隔采樣，不需要時間序列，只需要信號數組即可分析。

分析過程如下：

• 對于含有 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,4))
plt.plot(freqs,2*np.abs(xf),'r--')

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

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

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