目錄
- Python進行FFT頻譜分析
- FFT點數分析
- Cosine信号波形
- 周期方波信号波形
- 複合信号進行FFT(補零，截斷加長，加窗)
Python進行FFT頻譜分析
聲明：本文思想均來自陳愛軍老師《深入淺出通信原理》連載313-389
FFT點數分析
連載543
FFT點數 = OFDM符号周期 x 采樣頻率
OFDM符号周期 = 1/子載波間隔
Cosine信号波形
DFT公式: $X(k)=\frac{1}{N} \sum_{n=0}^{N-1}x(n)e^{-j\frac{2\pi}{N}kn} $
由常識可知\(cos(2\pi ft)\)的頻譜在1和-1處均有1/2的值，如下所示
即實數信号有對稱的頻譜形狀
# cosine wave(cos2pit) 進行 DFT
plt.subplot(2,1,1)
t = np.arange(0,1.01,0.01)
f = np.cos(2*np.pi*t)
plt.plot(t,f)
plt.xlabel('t')
plt.ylabel('magnitude')
plt.title('cosine wave')
N = 10
n = np.arange(0,N,1)
f = np.cos(2*np.pi*n/N)
t = np.arange(0,1,0.1)
plt.stem(t,f,use_line_collection=True)
# N=10, Ts=0.1s
plt.subplot(2,1,2)
n = np.arange(0,N,1)
x = np.cos(2*np.pi*n/N)
c = []
for k in np.arange(-N/2,N/2+1,1):
S=0
for n2 in n:
S = S+x[n2]*np.exp(-1j*2*np.pi*k*n2/N)
c.append(S/N)
k = np.arange(-N/2,N/2+1,1)
plt.stem(k,c,use_line_collection=True)
plt.xlabel('f')
plt.ylabel('Magnitude')
plt.title('DFT of cosine wave')
plt.tight_layout() 
直接調用FFT函數計算：
# 時域采樣信号
plt.subplot(211)
N = 10
n = np.arange(0,N,1)
f = np.cos(2*np.pi*n/N)
t = np.arange(0,1,0.1)
plt.stem(t,f,use_line_collection=True)
plt.xlabel('t')
plt.ylabel('magnitude')
plt.title('cosine wave')
# 頻譜信号
# REF: http://localhost:8888/notebooks/PycharmProjects/Untitled2.ipynb
plt.subplot(212)
S = np.fft.fft(f)/N
# 給出橫坐标的數字頻率, 第一個參數n是FFT的點數，一般取FFT之後的資料的長度（size）, 第二個參數d是采樣周期，其倒數就是采樣頻率Fs，即d=1/Fs
n = np.fft.fftfreq(N, 1/N)
# 用于将FFT變換之後的頻譜顯示範圍從[0, N]變為：[-N/2, N/2-1](N為偶數) 或者 [-(N-1)/2, (N-1)/2](N為奇數)
S = np.fft.fftshift(S)
n = np.fft.fftshift(n)
plt.stem(n,S,use_line_collection=True)
plt.xlabel('f')
plt.ylabel('Magnitude')
plt.title('DFT of cosine wave')
plt.tight_layout() 不限制頻域範圍[-2N,2N]，可見DFT後X(k)是周期為N的函數，是以通常取值區間取主值區間[-N/2,N/2]
# N=10, Ts=0.1s
plt.subplot(2,1,2)
n = np.arange(0,N,1)
x = np.cos(2*np.pi*n/N)
c = []
for k in np.arange(-N*2,N*2+1,1):
S=0
for n2 in n:
S = S+x[n2]*np.exp(-1j*2*np.pi*k*n2/N)
c.append(S/N)
k = np.arange(-N*2,N*2+1,1)
plt.stem(k,c,use_line_collection=True)
plt.xlabel('f')
plt.ylabel('Magnitude')
plt.title('DFT of cosine wave')
plt.tight_layout() IDFT變換的本質是将信号的樣本序列（N個樣點）表示成1個直流樣點序列，2/N個逆時針旋轉的旋轉向量樣點序列和2/N-1個順時針樣點序列
k=9的樣點資料與k=-1相同，因為\(e^{j\frac{2\pi}{N} kn}\)當 k=N/2+i 和 k=-N/2+i 相同
周期方波信号波形
# 周期方波作FFT
plt.figure(figsize=(6,10))
plt.subplot(511)
N = 40
x1 = np.ones(int(N/4))
x2 = np.append(x1,np.zeros(int(N/2)+1))
x3 = np.append(x2,np.ones(int(N/4)-1))
n = np.arange(0,len(x3),1)
plt.stem(n,x3,use_line_collection=True) # 偶對稱的周期方波
plt.title("偶對稱的周期方波",fontproperties='stsong')
plt.subplot(512)
X = np.fft.fft(x3)/N
plt.stem(n,X,use_line_collection=True) # 預設隻畫實部
plt.title("偶對稱的周期方波FFT實部",fontproperties='stsong')
# 虛部
plt.subplot(513)
plt.stem(n,np.imag(X),use_line_collection=True)
plt.title("偶對稱的周期方波FFT虛部",fontproperties='stsong')
# IFFT
# plt.subplot(5,1,3)
x = np.fft.ifft(X)
# plt.stem(n,x,use_line_collection=True)
# plt.title("IFFT",fontproperties='stsong')
# 幅度譜, 和連載347的不同
plt.subplot(514)
plt.stem(n,abs(X),use_line_collection=True)
plt.title("偶對稱的周期方波FFT幅度譜",fontproperties='stsong')
# 相位譜
plt.subplot(515)
plt.stem(n,np.angle(X),use_line_collection=True)
plt.title("偶對稱的周期方波FFT相位譜",fontproperties='stsong')
plt.tight_layout() 複合信号進行FFT(補零，截斷加長，加窗)
\(f(t) = cos(200 \pi t)+sin(100 \pi t)+cos(50 \pi t)\) 進行頻譜分析
# 連載354
fs = 400 #采樣頻率
T = 0.04 #截取總時間
N = fs*T #采樣16個點
t = np.arange(0,T,1/fs)
x = np.cos(200*np.pi*t)+np.sin(100*np.pi*t)+np.cos(50*np.pi*t)
plt.subplot(311)
plt.stem(t,x,use_line_collection=True)
plt.title('$f(t) = cos(200 \pi t)+sin(100 \pi t)+cos(50 \pi t)$')
plt.subplot(312)
X = np.fft.fft(x)
# n = np.arange(0,fs,fs/N) # 16個點，頻譜寬度400Hz
n = np.fft.fftfreq(int(N),1/fs)
X = np.fft.fftshift(X)
n = np.fft.fftshift(n)
plt.stem(n,abs(X),use_line_collection=True)
plt.title('Magnitude Spectrum')
plt.tight_layout() 顯然該信号頻譜很完美，但若要提高頻譜密度如何做？
答案：補零，将16個采樣點補零至256個采樣點
其實質是增大信号周期
# 采樣資料後補0來提高頻譜密度，即增大時域信号周期
fs = 400 #采樣頻率
T = 0.04 #截取總時間
N = fs*T #采樣16個點
t = np.arange(0,T,1/fs)
x = np.cos(200*np.pi*t)+np.sin(100*np.pi*t)+np.cos(50*np.pi*t)
t2 = np.arange(T,0.64,1/fs)
t = np.append(t,t2)
x = np.append(x,np.zeros(len(t2))) #補0增大頻譜密度
plt.subplot(311)
# plt.stem(t,x,use_line_collection=True)
plt.plot(t,x)
plt.title("補零增加頻譜密度的f(t)",fontproperties='stsong')
plt.subplot(312)
X = np.fft.fft(x)
n = np.arange(0,400,fs/N/N) # 16個點，頻譜寬度400Hz
plt.stem(n,abs(X),use_line_collection=True) # 顯然此時頻譜分辨率降低
plt.title("有旁瓣洩露的頻譜",fontproperties='stsong')
# plt.plot(n,abs(X))
plt.tight_layout() 頻譜密度提高了，但是頻譜分辨率顯然下降，出現旁瓣洩露
那麼如何提高分辨率呢？很明顯，增加截取長度即可
# 增加信号周期數提高頻譜分辨率
fs = 400 #采樣頻率
T = 0.16 #截取總時間，四個周期
N = fs*T #采樣64（Not 16）個點
t = np.arange(0,T,1/fs)
x = np.cos(200*np.pi*t)+np.sin(100*np.pi*t)+np.cos(50*np.pi*t)
t2 = np.arange(T,0.64,1/fs)
t = np.append(t,t2)
x = np.append(x,np.zeros(len(t2))) #補0增大頻譜密度
plt.subplot(311)
# plt.stem(t,x,use_line_collection=True)
plt.plot(t,x)
plt.title("增加4倍信号截取長度增加頻譜分辨率的f(t)",fontproperties='stsong')
plt.subplot(312)
X = np.fft.fft(x)
n = np.arange(0,fs,fs/256) # 32個點，頻譜寬度400Hz
# plt.stem(n,abs(X),use_line_collection=True) # 顯然此時頻譜分辨率降低
plt.plot(n,abs(X)[:256]) # 此時存在的非信号頻譜分量就是洩露效應
plt.title("分辨率提高的頻譜)",fontproperties='stsong')
plt.tight_layout() 除了原信号确實存在的頻率分量外還存在原信号中沒有的頻率分量，這就是洩露效應
洩露效應與截斷長度相關，截取長度越長，分辨率越高，洩露越少
截斷的過程就是與矩形窗相乘的過程，我們分析的頻譜是原始信号與矩形窗信号乘積再補零的頻譜
plt.figure(figsize=(6,6))
# 原信号f(t)
plt.subplot(3,1,1)
fs = 400 #采樣頻率
T = 0.04*8 #周期數目
N = fs*T #采樣點數
t = np.arange(-T/4,3*T/4,1/fs)
f = np.cos(200*np.pi*t)+np.sin(100*np.pi*t)+np.cos(50*np.pi*t)
plt.stem(t,f,use_line_collection=True)
plt.plot(t,f,color='r')
plt.title('$f(t) = cos(200 \pi t)+sin(100 \pi t)+cos(50 \pi t)$')
# 矩形窗w(t)
plt.subplot(312)
w1 = np.zeros(int(N/4))
w2 = np.ones(int(N/2))
w = np.append(w1,w2)
w = np.append(w,w1)
plt.stem(t,w,use_line_collection=True)
plt.title('Rectangle Windows w(t)')
# 乘積信号
plt.subplot(313)
y = f*w
plt.stem(t,y,use_line_collection=True)
plt.title('f(t)*w(t)')
plt.tight_layout() 連載367
截取時長等于頻率成分周期的整數倍，對采樣序列進行周期性擴充得到的周期序列的周期正好等于頻率成分的周期，可確定不出現洩露
采樣頻率也必須是信号頻率的整數倍，否則也會造成信号洩露
且如上所示補零也會造成頻譜洩露
FFT頻譜分析Python進行FFT頻譜分析
以對數角度來看，主瓣和旁瓣差别不大
plt.figure(figsize=(6,6))
# 原信号f(t)
plt.subplot(3,2,1)
fs = 400 #采樣頻率
T = 0.04*16 #周期數目
N = fs*T #采樣點數
t = np.arange(0,T,1/fs)
f = np.cos(200*np.pi*t)
t = np.arange(0,N,1)
plt.stem(t,f,use_line_collection=True)
# plt.plot(t,f,color='r')
plt.title('$f(t) = cos(200 \pi t)$')
# 矩形窗w(t)
plt.subplot(323)
w1 = np.ones(int(N/8))
w2 = np.zeros(int(N-N/8))
w = np.append(w1,w2)
plt.stem(t,w,use_line_collection=True)
plt.title('Rectangle Windows w(t)')
# 乘積信号
plt.subplot(325)
y = f*w
plt.stem(t,y,use_line_collection=True)
plt.title('y(t) = f(t)*w(t)')
# Spetrum
plt.subplot(322)
F = np.fft.fft(f)
plt.stem(t,F,use_line_collection=True)
plt.title('spectrum of f(t)')
plt.subplot(324)
W = np.fft.fft(w)
plt.stem(t,W,use_line_collection=True)
plt.title('Spectrum of w(t)')
plt.subplot(326)
Y = np.fft.fft(y)
# plt.stem(t,Y,use_line_collection=True)
plt.plot(t,Y)
plt.title('Spectrum of y(t)')
plt.tight_layout() plt.figure(figsize=(6,6))
# 原信号f(t)
plt.subplot(3,1,1)
fs = 400 #采樣頻率
T = 0.04*2 #周期數目
N = fs*T #采樣點數
t = np.arange(0,T,1/fs)
f = np.cos(200*np.pi*t)
t2 = np.arange(T,0.64,1/fs)
t = np.append(t,t2)
f = np.append(f,np.zeros(len(t2))) #補0增大頻譜密度
plt.plot(t,f)
plt.title('$f(t) = cos(200 \pi t))$')
plt.subplot(312)
F = np.fft.fft(f)
plt.plot(t,F)
plt.title('Spectrum')
plt.subplot(313)
plt.semilogy(t,abs(F))
# Z = 20*np.log10(abs(F))
# plt.plot(t,Z)
plt.ylim([0.001,100])
plt.title('Spectrum in semilogy')
plt.xlabel('Frequency(Hz)')
plt.ylabel('Aplitude(dB)')
plt.tight_layout() \(w_hnn(n）=0.5-0.5cos(\frac{2\pi n}{N-1}), n=0, 1, 2, \cdots, N-1\)
漢甯窗主瓣寬度增大，旁瓣高度減小；即犧牲分辨率，換來旁瓣洩露的減小
# hanning window
plt.subplot(211)
N = 256
y1 = np.hanning(32)
y2 = np.zeros(int(N-32))
y = np.append(y1,y2)
n = np.arange(0,N,1)
plt.stem(n,y,use_line_collection=True)
plt.title('hanning window')
# Spectrum
plt.subplot(212)
Y = np.fft.fft(y)
# plt.stem(n,abs(Y)/abs(len(Y)),use_line_collection=True)
# plt.plot(n,abs(Y)/max(abs(Y)))
plt.semilogy(n,abs(Y)/max(abs(Y)))
plt.tight_layout() \(w_{hmm}(n)=0.54-0.46cos(\frac{2\pi n}{N-1}), n=0, 1, 2, \cdots, N-1\)
漢明窗相對矩形窗主瓣寬度增大，旁瓣高度減小；即犧牲分辨率，換來旁瓣洩露的減小。相對漢甯窗旁瓣高度減小不明顯，衰減到最大幅度的8%左右而不是衰減至零，Blackman窗則主瓣更大，旁邊更小。
# 窗函數對比
N = 4096
n = 32
y1a = np.ones(n)
y1b = np.hanning(n)
y1c = np.hamming(n)
y1d = np.blackman(n)
y2 = np.zeros(N-n)
n = np.arange(0,N,1)
# d = {'a':'rectangle', 'b':'hanning', 'c':'hamming', 'd':'blackman'}
d = ['rectangle', 'hanning', 'hamming', 'blackman']
count = 0
for y in [y1a,y1b,y1c,y1d]:
yx = np.append(y,y2)
Yx = np.fft.fft(yx)
plt.semilogy(n,abs(Yx)/max(abs(Yx)),label=d[count])
count+=1
plt.ylim(0.000001,1)
plt.xlim(0,N/2)
plt.legend(loc='upper right')
plt.yticks([1,0.1,0.01,0.001,0.0001,0.00001,0.000001],['0','-20','-40','-60','-80','-100','-120'])
plt.ylabel('Aplitude(dB)')
plt.title('窗函數比較',fontproperties='stsong') plt.figure(figsize=(12,6))
# 原信号f(t)
plt.subplot(231)
fs = 400 #采樣頻率
T = 0.04*4 #周期數目
N = 256 #采樣點數
t = np.arange(0,T,1/fs)
f = np.cos(200*np.pi*t)+np.sin(100*np.pi*t)+np.cos(50*np.pi*t)
N0 = int(N - T*fs)
ft = np.append(f,np.zeros(N0))
t = np.arange(0,fs,fs/N)
plt.plot(t,ft,color='r')
plt.title('$f(t) = cos(200 \pi t)+sin(100 \pi t)+cos(50 \pi t)$')
# add rectangle window
plt.subplot(234)
X = np.fft.fft(ft)
plt.plot(t,abs(X),label='rectangle')
plt.title('Retangle window Spectrum')
# add hamming window
plt.subplot(232)
x1c = np.hamming(T*fs)
ft = f*x1c
ft = np.append(ft,np.zeros(N0))
plt.plot(t,ft)
plt.title('$f(t)=f(t)*hamming(T*fs)$')
plt.subplot(235)
Xc = np.fft.fft(ft)
plt.plot(t,Xc)
plt.title('Hamming Window Spectrum')
# add blackman window
plt.subplot(233)
x1d = np.blackman(T*fs)
ft = f*x1d
ft = np.append(ft,np.zeros(N0))
plt.plot(t,ft)
plt.title('$f(t)=f(t)*blackman(T*fs)$')
plt.subplot(236)
Xd = np.fft.fft(ft)
plt.plot(t,Xd)
plt.title('blackman Window Spectrum')
plt.tight_layout() 照慣例附上Github倉庫對應代碼連結
轉載于:https://www.cnblogs.com/WindyZ/p/FFT_Analysis.html