DSP筆記1_使用 NumPy自帶的fft做_觀察fft結果、ifft後的訊號波形圖、頻率對應相位
在Python當中我們可使用 FFT 函式，將時域訊號轉換成頻域訊號。
例如使用 NumPy自帶的：
import numpy.fft as fft
而另一套scipy引入方式則是用
import scipy.fftpack as scifft
- fft功能:主要就是把時間序列訊號做FFT，獲得其幾個頻率、振福和初始相位。
- ifft功能:也就是把複數值序列做傅立葉逆變化，獲得複數值的訊號，取其實部，是時頻訊號。可以寫出訊號的余弦函數的結構式。
以下就DSP幾個術語做個介紹
訊號頻率：
指訊號本身週期性變化的速度，通常使用 f 表示，單位為赫茲（Hz）。
例如，一秒鐘重複 10 次的訊號，其頻率為 10 Hz。
取樣頻率：
指每秒鐘對連續訊號取樣的次數，通常使用 Fs 表示，單位也是 Hz。
表示每秒取得 100 個取樣點。
奈奎斯(Nyquist-Shannon sampling theorem)取樣定理：
為了正確還原訊號，取樣頻率至少必須大於訊號最高頻率的兩倍。
舉例:像是人類聽力範圍坐落於20赫茲-20000赫茲(20kHz)
根據奈奎斯取樣定理，fs > 2 fH 也就是 fs > 2 * 20kHz
目前網際網路流傳的音樂檔mp3，其取樣頻率也通常是設在44kHz、48kHz。
實際應用通常會選擇比兩倍更高的取樣頻率，以降低混疊現象。
FFT 頻譜幅值：
fft_magnitude = np.abs(np.fft.fft(signal, n=NFFT))
np.abs() 用來計算 FFT 複數結果的絕對值，也就是頻譜的幅值。
振幅(頻)譜(amplitude spectrum):
透過傅立葉變換（Fourier Transform）將時域訊號轉成頻域後，取其絕對值或強度大小而得到的結果。
能量訊號：
若訊號振幅的絕對值平方，在整個時間範圍 (−∞,+∞) 內的積分為有限值，則稱為能量訊號。
常見例子包括有限時間的方波訊號、三角形訊號、脈衝訊號、暫態訊號、非週期性的確定訊號，以及非隨機訊號。
能量頻譜：
又稱為「能量頻譜密度」。圖形的橫軸表示頻率，縱軸表示訊號在各頻率上的能量。
訊號總能量可由各頻率成分的幅值平方加總或積分得到。
實際分析時，可以利用快速傅立葉轉換 FFT 計算訊號的頻率成分，再將頻譜幅值取絕對值平方，得到能量頻譜。
能量頻譜的實際意義：
能量頻譜用來表示訊號總能量分布於不同頻率上的情形。透過能量頻譜，可以觀察訊號的主要頻率成分，以及各頻率成分所包含的能量大小。
功率訊號：
若訊號的絕對值平方，在時間區間 (−T/2，T/2) 內的積分平均值，於 T 趨近無限大時存在有限極限，則稱為功率訊號。
常見例子包括週期訊號、常數訊號、階躍訊號及隨機訊號。
功率頻譜、譜密度(PSD,power spectral density)：
又稱為「功率頻譜密度」。圖形的橫軸表示頻率，縱軸表示功率，用來描述訊號在單位頻寬(頻帶)內所具有的功率。
在時域中，可將訊號各取樣值的絕對值平方加總，再除以取樣點總數，以估算訊號的平均功率。
在實際工程應用中，不能只直接使用 FFT 結果；通常需要將 FFT 的絕對值平方，
再依訊號持續時間 T、取樣點數或取樣頻率進行正規化，才能得到具有正確尺度的功率頻譜或功率頻譜密度。
功率頻譜的實際意義：
功率頻譜用來表示訊號功率分布於各個頻率上的情形。
藉由功率頻譜，可以觀察訊號的主要頻率成分，以及各頻率成分所占的功率大小。
以下用這個正弦波來去做模擬訊號的fft頻譜分析
練習程式(使用 NumPy自帶的fft)
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 | import matplotlib.pyplot as plt import numpy as np import numpy.fft as fft plt.rcParams['font.family'] = 'sans-serif' plt.rcParams['font.sans-serif'] = [ 'Microsoft JhengHei', # 微軟正黑體 'Noto Sans CJK TC', 'Arial Unicode MS' ] plt.rcParams['axes.unicode_minus'] = False Fs = 256 #採樣頻率：256 Hz，1秒採樣256個點。 dt = 1/Fs #採樣間隔、採樣週期(時間間隔、步長)，也是採樣頻率的倒數。兩個採樣點之間時間間隔。 L = 1024 #訊號長度(採樣點數)，使用2的次方 2^10 t = [i*dt for i in range(L)] #從0開始到L-1個數值 t = np.array(t)#把list轉化給numpy中的array #偏移量=2、振幅 = 3、頻率 = 80 Hz、初始相位 = -60度 S = 2+3*np.sin(2*np.pi*80*t-60*np.pi/180)#初始相位(角度)、振幅、頻率 complex_array = fft.fft(S) print(complex_array) print(complex_array.shape)#(1024,) print(complex_array.dtype)#complex128 print(complex_array[1]) N_Upper=100 #========================5-1========================== #原始訊號跟它的fft圖像、ifft圖像 plt.figure(figsize=(16,8))#寬:16,高:8 plt.subplot(511)#繪製5橫列、一直行 plt.grid(linestyle=':') #plt.plot(橫軸時間變量0~49前50個,縱軸訊號變量0~49前50個,label='S(t)') plt.plot(t[0:N_Upper],S[0:N_Upper],label='S(t)') plt.xlabel('t(秒)') plt.ylabel('S(t)') plt.title('原始訊號圖') plt.legend() #========================5-2========================== #原始訊號的fft實數部分圖 plt.subplot(512)#繪製5橫列、一直行 plt.grid(linestyle=':') plt.plot(t[0:N_Upper],complex_array.real[0:N_Upper],label='fft實部') plt.xlabel('t(秒)') plt.ylabel('fft實部') plt.title('原始訊號fft實部圖') plt.legend() #========================5-3========================== #原始訊號的fft虛數部分圖 plt.subplot(513)#繪製5橫列、一直行 plt.grid(linestyle=':') plt.plot(t[0:N_Upper],complex_array.imag[0:N_Upper],label='fft虛部') plt.xlabel('t(秒)') plt.ylabel('fft虛部') plt.title('原始訊號fft虛部圖') plt.legend() #========================5-4========================== #原始訊號fft之後結果再做ifft變化序列 plt.subplot(514)#繪製5橫列、一直行 ifft_s = fft.ifft(complex_array)#ifft_s會是複數 plt.grid(linestyle=':') plt.plot(t[0:N_Upper],ifft_s.real[0:N_Upper],label='ifft_s',color='orangered') plt.xlabel('t(秒)') plt.ylabel('ifft_s(t)幅值') plt.title('ifft轉換圖') plt.legend() #========================5-5========================== #獲得分解波的頻率序列 #t.size ->時間變量的長度 #t[1]-t[0] ->採樣間隔(採樣率倒數)，相鄰採樣點時間距離。 freqs = fft.fftfreq(t.size,t[1]-t[0]) #FFT複數結果(Magnitude)幅度:實部的平方+虛部的平方再去開根號 magni = np.abs(complex_array)#FFT複數的模長(未正規化的雙邊 FFT模長) plt.subplot(515)#繪製5橫列、一直行 plt.xlabel('頻率(Frequency)') plt.ylabel('FFT幅度值') plt.title('FFT變化(頻幅度圖)') plt.tick_params(labelsize=10)#坐標軸標籤文字大小 plt.grid(linestyle=':') #繪製正頻率圖像->跟原始訊號振幅=3有落差，已經超過1000了。 plt.plot(freqs[freqs > 0],magni[freqs > 0],c='orangered',label='Frequency') plt.legend() plt.tight_layout()#避免圖表標籤有重疊 plt.show() #========================1-1========================== plt.figure(figsize=(16,8))#寬:16,高:8 angle_y = np.arctan2(complex_array.imag,complex_array.real) plt.plot(freqs[freqs>0],angle_y[freqs>0],c='orangered',label='Frequency') plt.xlabel('頻率(Frequency)') plt.ylabel('相位') plt.legend() plt.tick_params(labelsize=10)#坐標軸標籤文字大小 plt.show() |
運行後第一張圖呈現的是五個chart合在同一個圖表呈現
可觀察到fft結果經過ifft後的訊號波形圖幾乎吻合。
運行後第二張圖呈現的是頻率對應相位
留言
張貼留言