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

可是自己動手實現(xiàn)一遍才是最好的學習。

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

分析過程如下：

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

使用方法：

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

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

完整程序：

# 離散時間傅里葉變換的 python 實現(xiàn)
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()

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

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

最新評論