提高 FFT 分辨率 - 缩放 Y 轴
Improving FFT resolution - scaling Y axis
我正在尝试绘制 Python 中 FM 信号的贝塞尔等效值的频域。
我遇到的第一个问题是提高 FFT 的分辨率以查看更小的带宽 - 我相信我已经通过将贝塞尔波函数的总和与一个大小为保存波形数据的数组的两倍。
但是,这对 Y 轴的缩放有影响 - 当我增加零填充的大小时,Y 轴的大小会下降。
我可以将 Y 轴乘以什么因子来反转它?到目前为止的实验让我走到了死胡同...
非常感谢!
import numpy as np
import matplotlib.pyplot as plt
import math as mt
import scipy.special as sc
def bWave(fc=10000, B=1.25):
time = np.linspace(0, 0.5, 20001, True)
data = np.zeros(len(time))
retarr = np.column_stack((time,data))
Vc = 6
fm = 3
for n in range(-5,5):
for row in range(len(retarr)):
retarr[row][1] += Vc*sc.jv(n,B)*mt.cos(2*mt.pi*
(fc+n*fm)*retarr[row][0])
return retarr
FM_array = bWave()
# ------------- SIGNAL PLOT OF AM -----------------
scaling = 2 #default 2, cuts out symmetry from FFT
buffer_ratio = 1
padded_array =
np.concatenate((FM_array[:,1],np.zeros(buffer_ratio*len(FM_array[:,1])))) #pad array with zeros
Y = np.fft.fft(padded_array) #perform FFT
N = len(Y)/scaling + 1 # get FFT length (avoid reflection)
T = FM_array[1][0] - FM_array[0][0] #get time interval of FFT
fa = 1/T #sampling frequency
Xaxis = np.linspace(0, fa/scaling, N, endpoint=True) # create x axis vector from 0 to nyquist freq. (fa/2) with N values
plt.plot(Xaxis, (2.0/((N)*scaling)) * np.abs(Y[0:N])) # multiply y axis by 2/N to get actual values
plt.grid(True)
plt.show()
几点:
- 您确定
bWave()
函数正确吗?您的贝塞尔函数不依赖于时间,因此您可以通过傅里叶变换余弦轻松找到封闭形式的解决方案。
- 不要用零填充,而是增加
bWave()
信号的时间周期(见下面的代码)以提高频率分辨率。
- 使用
numpy
而不是 math
函数。它使您的代码更具可读性和更快。
以下代码绘制了不同时间段的 FFT。时间周期越长,峰值就会变得越尖锐(余弦的傅里叶变换):
import numpy as np
import matplotlib.pyplot as plt
import scipy.special as sc
def bWave2(t, fc=10000, B=1.25):
""" Useing only numpy """
Vc, fm = 6, 3
y = np.zeros(len(t))
for n in range(-5,5):
y += Vc*sc.jv(n,B)*np.cos(2*np.pi*(fc+n*fm)*t)
return y
fg, ax = plt.subplots(1, 1)
fc=10000
for q in range(0, 5):
k = 15001*(1+q)
t = np.linspace(0-0.25*q, 0.5+0.25*q, k)
y = bWave2(t, fc)
Y = np.fft.rfft(y)/N # since y is real, rfft() instead of fft is used
f = np.fft.rfftfreq(k, t[1] - t[0]) # frequencies for rfft()
ax.plot(f, np.abs(Y), label=u"$\tau=${}".format(t[-1]-t[0]), alpha=.5)
ax.set_xlim(fc-50, fc+50)
ax.grid(True)
ax.legend(loc="best")
fg.canvas.draw()
plt.show()
请注意 rfft(y)/N
中的 /N
只是为了具有可比较的 FFT 值而添加的。由于采样率是恒定的,能量,即 |Y(f)|²,随着时间的增加而增加。在您的代码中,采样率发生了变化,信号的能量也发生了变化。
我正在尝试绘制 Python 中 FM 信号的贝塞尔等效值的频域。
我遇到的第一个问题是提高 FFT 的分辨率以查看更小的带宽 - 我相信我已经通过将贝塞尔波函数的总和与一个大小为保存波形数据的数组的两倍。
但是,这对 Y 轴的缩放有影响 - 当我增加零填充的大小时,Y 轴的大小会下降。 我可以将 Y 轴乘以什么因子来反转它?到目前为止的实验让我走到了死胡同...
非常感谢!
import numpy as np
import matplotlib.pyplot as plt
import math as mt
import scipy.special as sc
def bWave(fc=10000, B=1.25):
time = np.linspace(0, 0.5, 20001, True)
data = np.zeros(len(time))
retarr = np.column_stack((time,data))
Vc = 6
fm = 3
for n in range(-5,5):
for row in range(len(retarr)):
retarr[row][1] += Vc*sc.jv(n,B)*mt.cos(2*mt.pi*
(fc+n*fm)*retarr[row][0])
return retarr
FM_array = bWave()
# ------------- SIGNAL PLOT OF AM -----------------
scaling = 2 #default 2, cuts out symmetry from FFT
buffer_ratio = 1
padded_array =
np.concatenate((FM_array[:,1],np.zeros(buffer_ratio*len(FM_array[:,1])))) #pad array with zeros
Y = np.fft.fft(padded_array) #perform FFT
N = len(Y)/scaling + 1 # get FFT length (avoid reflection)
T = FM_array[1][0] - FM_array[0][0] #get time interval of FFT
fa = 1/T #sampling frequency
Xaxis = np.linspace(0, fa/scaling, N, endpoint=True) # create x axis vector from 0 to nyquist freq. (fa/2) with N values
plt.plot(Xaxis, (2.0/((N)*scaling)) * np.abs(Y[0:N])) # multiply y axis by 2/N to get actual values
plt.grid(True)
plt.show()
几点:
- 您确定
bWave()
函数正确吗?您的贝塞尔函数不依赖于时间,因此您可以通过傅里叶变换余弦轻松找到封闭形式的解决方案。 - 不要用零填充,而是增加
bWave()
信号的时间周期(见下面的代码)以提高频率分辨率。 - 使用
numpy
而不是math
函数。它使您的代码更具可读性和更快。
以下代码绘制了不同时间段的 FFT。时间周期越长,峰值就会变得越尖锐(余弦的傅里叶变换):
import numpy as np
import matplotlib.pyplot as plt
import scipy.special as sc
def bWave2(t, fc=10000, B=1.25):
""" Useing only numpy """
Vc, fm = 6, 3
y = np.zeros(len(t))
for n in range(-5,5):
y += Vc*sc.jv(n,B)*np.cos(2*np.pi*(fc+n*fm)*t)
return y
fg, ax = plt.subplots(1, 1)
fc=10000
for q in range(0, 5):
k = 15001*(1+q)
t = np.linspace(0-0.25*q, 0.5+0.25*q, k)
y = bWave2(t, fc)
Y = np.fft.rfft(y)/N # since y is real, rfft() instead of fft is used
f = np.fft.rfftfreq(k, t[1] - t[0]) # frequencies for rfft()
ax.plot(f, np.abs(Y), label=u"$\tau=${}".format(t[-1]-t[0]), alpha=.5)
ax.set_xlim(fc-50, fc+50)
ax.grid(True)
ax.legend(loc="best")
fg.canvas.draw()
plt.show()
请注意 rfft(y)/N
中的 /N
只是为了具有可比较的 FFT 值而添加的。由于采样率是恒定的,能量,即 |Y(f)|²,随着时间的增加而增加。在您的代码中,采样率发生了变化,信号的能量也发生了变化。