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

可是自己動(dòng)手實(shí)現(xiàn)一遍才是最好的學(xué)習(xí)。

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

分析過(guò)程如下：

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

使用方法：

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

• 計(jì)算出的每一個(gè)頻率下的幅值需要除以時(shí)間序列的長(zhǎng)度，類似求平均的過(guò)程
• 每一個(gè)頻率下的幅值是一個(gè)復(fù)數(shù)，需要對(duì)它求模，而且因?yàn)樵谪?fù)頻率處也有值，所以需要對(duì)于實(shí)信號(hào)需要乘 2
• 頻率的序列為 0 到采樣率的一半，長(zhǎng)度為 n/2+1
  

完整程序：

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

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

以上就是本文的全部?jī)?nèi)容，希望對(duì)大家的學(xué)習(xí)有所幫助，也希望大家多多支持腳本之家。

最新評(píng)論