我有一个以 2MHz 采样率记录的信号。在可能需要任何抽取之前，我首先使用 STFT/频谱图在我的周期性记录信号中寻找峰值。从这里，我可以清楚地看到一些山峰。使用 PSD，我发现在 ~500kHz 之后频率并不突出，因此我应用巴特沃斯滤波器滤除高于此的频率。此外，我的天线从 ~20kHz 开始记录，因此低于此的频率也会被过滤掉。
该sample.wav文件可以在这里下载：https ://www.dropbox.com/s/hcwrtuzqbipatri/sample.wav?dl=0
滤波器幅度响应如下所示：
然而，比较原始信号和带通滤波信号的频谱图看起来完全不同，我似乎无法理解为什么。原始信号在左边，巴特沃斯滤波信号在右边。y 轴刻度是频率 (kHz)。
过滤器使用的代码如下：
import numpy as np
import os
import scipy.io.wavfile
from scipy import signal
# import librosa
import matplotlib.pyplot as plt
from sklearn.decomposition import PCA
fs, x = scipy.io.wavfile.read('sample.wav') # fs is 2MHz (2000000)
# Convert to mono by avg I (left) and Q (right) channels
x = np.mean(x, axis=1)
x = np.asarray(x, dtype='float')
# Once the recording is in memory, we normalise it to +1/-1
# x /= np.max(np.abs(x)) ## This is sensitive to outliers and rescaling is not consistent
x /= x.std()
t_tot = 1 / fs*x.shape[0]
t = np.linspace(0, t_tot, x.shape[0])
tmin = np.min(t)
tmax = np.max(t)
N = len(x)
f_l, f_u = 20000, 500000 # Band from 20kHz to 500kHz
wp = np.array([f_l, f_u])*2/fs # normalized pass band frequnecies
ws = np.array([0.8*f_l, 1.2*f_u])*2/fs # normalized stop band frequencies
b, a = signal.iirdesign(wp, ws, gpass=60, gstop=80, ftype="butter", analog=False)
# plot filter response:
w, h = signal.freqz(b, a, whole=False)
ff_w = w*fs/(2*np.pi)
fg, ax = plt.subplots()
ax.set_title('Butterworth filter amplitude response')
ax.plot(ff_w/1000, np.abs(h))
ax.set_ylabel('Relative Amplitude')
ax.grid(True)
ax.set_xlabel('Frequency (kHz)')
fg.canvas.draw()
# do the filtering:
zi = signal.lfilter_zi(b, a)*x[0]
x1, _ = signal.lfilter(b, a, x, zi=zi)
# calculate the avarage:
avg = np.mean(x1**2)
print("RMS values is %g" % avg)
plt.show()
fig, ax = plt.subplots(1, 2, sharex=True, sharey=True)
f, t, Sxx = signal.spectrogram(x, fs=fs, window='hann', nperseg=8192, scaling='spectrum')
ax[0].pcolormesh(t, f / 1000, 10 * np.log10(Sxx), shading='gouraud')
f_1, t_1, Sxx_1 = signal.spectrogram(x1, fs=fs, window='hann', nperseg=8192, scaling='spectrum')
ax[1].pcolormesh(t_1, f_1 / 1000, 10 * np.log10(Sxx_1), shading='gouraud')
# plt.set_title(f'nperseg={str(i)}')
plt.xlabel('Time (s)')
plt.ylabel('')
plt.show()
