数据及问题描述:
我有一个以 1000 Hz 采样的信号。我在 120 Hz 处对其进行低通滤波,并希望制作低于此阈值的频率的频谱图。我正在使用 scipy 函数fftpack.fftfreq 来获取傅立叶系数,然后fftpack.fft 进行实际变换。
信号很长,我希望使用 50 毫秒窗口的频谱图长度约为 5 秒。我还过滤掉负频率进行绘图。给定此窗口中的样本数(50),我得到 50 个傅立叶系数。但是,其中一半是负数,另一部分介于 120 Hz 和 500 Hz 之间(自然,因为采样率为 1000 Hz)。
这给我留下了相当低分辨率的频谱图。我在感兴趣的频率区域(0、20、40、60、100 和 120 Hz)中只有大约 6 个块,然后有 18 个块在高达 500 Hz 的频率中显示低活动,因为它们被过滤掉了。
问题:
例如,我怎么能在 FFT 中指定 0 到 120 Hz 之间的 50 个频率区间?
我尝试了什么:
np.linspace我尝试通过使用来指定频率(而不是)来做这样的事情fft.fftfreq,但这引入了一些其他错误,即频谱图总是看起来像一面镜子,在最高和最低频率处具有高功率,在中间具有低功率图表,与范围无关。老实说,我不确定为什么 linspace 而不是 fftfreq 会发生这种情况。它们都返回浮点数组。
代码示例包括在下面。提前感谢您的帮助!干杯
X = samples_filt # long 1-D vector of low-pass filtered samples
fs = 1000
time = .05 # window length in sec
N = int(fs * time) # num samples
tot_len = 5 * fs # 5 sec of the whole signal
X = X[:tot_len]
f = fftpack.fftfreq(N, 1.0/fs)
# f = np.ceil(np.linspace(0, max_freq, 52)[1:51]) # introduced error described above
mask = (f > 0) # mask for positive freqs
n_max = int(np.ceil(X.shape[0] / N)) # the number of segments of length N in the sample array data
f_values = np.sum(1 * mask) # how many values meet mask reqs
spectogram_data = np.zeros((f_values, n_max))
window = sp.signal.blackman(N) # taper used to improve contrast of spectrogram
for n in range(0, n_max):
subdata = X[(N * n):(N * (n + 1))]
F = fftpack.fft(subdata * window)
spectogram_data[:, n] = np.log(abs(F[mask]))
fig, ax = plt.subplots(1, 1, figsize=(8, 6))
p = ax.imshow(spectogram_data, origin='lower',
extent=(0, X.shape[0] / fs, 0, max(f)),
aspect='auto',
cmap=mpl.cm.RdBu_r)
cb = fig.colorbar(p, ax=ax)
样本输出:
例如,如果有相同的图表,其中包含 24 个频率区间,最高可达 120 Hz,那就太好了。
