FFT頻譜分析

Python進行FFT頻譜分析

聲明：本文思想均來自陳愛軍老師《深入淺出通信原理》連載313-389

目錄

• Python進行FFT頻譜分析
  • FFT點數分析
  • Cosine信號波形
  • 周期方波信號波形
  • 復合信號進行FFT(補零，截斷加長，加窗)

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

截取時長等於頻率成分周期的整數倍，對采樣序列進行周期性擴展得到的周期序列的周期正好等於頻率成分的周期，可確保不出現泄露

采樣頻率也必須是信號頻率的整數倍，否則也會造成信號泄露

且如上所示補零也會造成頻譜泄露

以對數角度來看，主瓣和旁瓣差別不大

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倉庫對應代碼鏈接