數字訊號處理 觀察訊號的頻譜

今天學習將時域訊號通過fft轉換為頻域訊號之後，將其各個頻率分量的幅值繪製成圖，可以很直觀地觀察訊號的頻譜。重點理解fft變換的過程。

程式來自參考書《python科學計算》

import numpy as np

import pylab as pl
from pylab import mpl
#首先定義兩個常數：sampling_rate, fft_size，
sampling_rate =
8000
#取樣頻率
fft_size =
512#fft的長度
#呼叫np.arange產生1秒鐘的取樣時間，t中的每個數值直接表示取樣點的#時間，因此其間隔為取樣週期1/sampline_rate
t = np.arange(0,
1.0,
1.0/sampling_rate)
#兩個正弦波的疊加,注意此兩波頻率正好為計算波形週期的10倍和15倍
x = np.sin(
2*np.pi*
156.25
*t)+
2*np.sin(
2*np.pi*
234.375
*t)#從波形資料x中擷取fft_size個點進行fft計算
xs = x[
:fft_size]
#計算波形能量
xf = np.fft.rfft(xs)
/fft_size
#計算出返回值中每個下標對應的真正的頻率。(源程式報錯：將/改為//，取整)
freqs = np.linspace(
0, sampling_rate//
2, fft_size//2+
1,dtype=
'int'
)#防止0幅值的成分造成log10無法計算，呼叫np.clip對xf的幅值進行上下限處理
xfp =
20*np.log10(np.clip(np.
abs(xf),1e
-20,1e100))
#時域波形和頻域波形繪製
mpl.rcparams[
'font.sans-serif']=
['simhei'
]#microsoft yahei
mpl.rcparams[
'axes.unicode_minus']=
false
pl.figure(figsize=(8
,4))
pl.subplot(
211)
pl.plot(t[
:fft_size]
, xs)
pl.xlabel(u"時間(秒)"
)pl.title(u"156.25hz和234.375hz的波形和頻譜"
)pl.subplot(
212)
pl.plot(freqs, xfp)
pl.xlabel(u"頻率(hz)"
)pl.subplots_adjust(hspace=
0.4)
pl.show(
)

頻譜圖中兩峰值對應的幅值為：

>>

> xfp[10]
-6.0205999132796251
>>
> xfp[15]
-9.6432746655328714e-16

pl.plot(freqs, np.

abs(xf[
:257])
)

為什麼選擇這兩個奇怪的頻率呢？因為這兩個頻率的正弦波在512個取樣點中正好有整數個週期。滿足這個條件波形的fft結果能夠精確地反映其頻譜。假設取樣頻率為fs,

取波形中的n個資料進行fft變換。那麼這n點資料報含整數個週期的波形時，fft所計算的結果是精確的。於是能精確計算的波形的週期是:

n*fs/n。對於8khz取樣，512點fft來說，8000/512.0 =

15.625hz，前面的156.25hz和234.375hz正好是其10倍和15倍。

即156.25hz的成分為-6db， 而234.375hz的成分為0db，與波形的計算公式中的各個分量的能量(振幅值/2)符合。可以這樣理解，這兩頻率對應的幅值應該為對應的兩正弦波的幅值，即對應為1和2，而20log10(1)和20log10(2)得到的結果正好是0db和6db，而上面的程式只取了頻譜的一半進行能量計算，所以對應的結果為20log10(0.5)和20log10(1)，即-6db和0db。

將波形計算公式修改為：

x = np.sin(

2*np.pi*
200*t)+2
*np.sin(
2*np.pi*
300*t)

這次得到的頻譜不再是兩個完美的峰值，而是兩個峰值頻率周圍的頻率都有能量。這顯然和兩個正弦波的疊加波形的頻譜有區別。本來應該屬於200hz和300hz的能量分散到了周圍的頻率中，這個現象被稱為頻譜洩漏。出現頻譜洩漏的原因在於fft_size個取樣點無法放下整數個200hz和300hz的波形。

頻譜洩漏的解釋

我們只能在有限的時間段中對訊號進行測量，無法知道在測量範圍之外的訊號是怎樣的。因此只能對測量範圍之外的訊號進行假設。而傅利葉變換的假設很簡單：測量範圍之外的訊號是所測量到的訊號的重複。

現在考慮512點fft，從訊號中取出的512個資料就是fft的測量範圍，它計算的是這512個資料一直重複的波形的頻譜。顯然如果512個資料報含整數個週期的話，那麼得到的結果就是原始訊號的頻譜，而如果不是整數週期的話，得到的頻譜就是如下波形的頻譜，這裡假設對50hz的正弦波進行512點fft：

由於這個波形的前後不是連續的，出現波形跳變，而跳變處的有著非常廣泛的頻譜，因此fft的結果**現頻譜洩漏。

接下來就要引入窗函式了。

數字訊號處理（1） 頻譜分析

離散傅利葉變換 dft 是傅利葉變換在時域和頻域上都呈現離散的形式，將時域訊號的取樣變換為在頻域的取樣。在實際應用中通常採用快速傅利葉變換 fft 以高效計算dft。其中n為dft的點數，點數越大頻率解析度越高，k 0,1,2,n 1。離散傅利葉變換可以看做是離散時域訊號與不同頻率的離散正弦訊號進行...

數字訊號處理

1.乙個切比雪夫i型模擬帶阻濾波器用下面的指標設計 通帶截止頻率為 和 阻帶截止頻率為 和 峰值通帶紋波是 最小阻帶衰減為 相應的模擬低通濾波器的頻帶截止頻率和階次是多少？帶阻濾波器的階次是多少？用matlab函式cheblord來驗證濾波器階次結果 寫出程式關鍵步驟 matlab驗證 n 3 fs...

數字訊號處理

一 訊號處理的典型過程 1 模擬濾波 x t xa t 觀測訊號經過前置模擬器ha s 去掉一些帶外成分和干擾 2 取樣 xa t xa nt 以取樣週期t對xa t 進行取樣，得到時域離散訊號xa nt 將時間離散化 3 a d 變換 xa nt x n 把原先訊號幅值連續變換 量化幅值 將幅值離...