numpy - Scipy Butter 带通没有产生预期的结果
问题描述
所以我正在尝试对 wav PCM 24 位 44.1khz 文件进行带通滤波。我想做的是带通0Hz-22Khz的每个频率。
到目前为止,我已经加载了数据并可以在 Matplot 上显示它,如下所示。
但是当我去应用我从这里得到的带通滤波器时
http://scipy-cookbook.readthedocs.io/items/ButterworthBandpass.html
所以我试图通过 100-101Hz 的带通作为测试,这是我的代码:
from WaveData import WaveData
import matplotlib.pyplot as plt
from scipy.signal import butter, lfilter, freqz
from scipy.io.wavfile import read
import numpy as np
from WaveData import WaveData
class Filter:
def __init__(self, wav):
self.waveData = WaveData(wav)
def butter_bandpass(self, lowcut, highcut, fs, order=5):
nyq = 0.5 * fs
low = lowcut / nyq
high = highcut / nyq
b, a = butter(order, [low, high], btype='band')
return b, a
def butter_bandpass_filter(self, data, lowcut, highcut, fs, order):
b, a = self.butter_bandpass(lowcut, highcut, fs, order=order)
y = lfilter(b, a, data)
return y
def getFilteredSignal(self, freq):
return self.butter_bandpass_filter(data=self.waveData.file['Data'], lowcut=100, highcut=101, fs=44100, order=3)
def getUnprocessedData(self):
return self.waveData.file['Data']
def plot(self, signalA, signalB=None):
plt.plot(signalA)
if signalB != None:
plt.plot(signalB)
plt.show()
if __name__ == "__main__":
# file = WaveData("kick.wav")
# fileA = read("kick0.wav")
f = Filter("kick.wav")
a, b = f. butter_bandpass(lowcut=100, highcut=101, fs=44100)
w, h = freqz(b, a, worN=22000) ##Filted signal is not working?
f.plot(h, w)
print("break")
我不明白我哪里出错了。
谢谢
解决方案
@WoodyDev 说的是真的:44.1 kHz 中的 1 Hz对于任何类型的滤波器来说都太小了。 只需查看滤波器系数返回:butter
In [3]: butter(5, [100/(44.1e3/2), 101/(44.1e3/2)], btype='band')
Out[3]:
(array([ 1.83424060e-21, 0.00000000e+00, -9.17120299e-21, 0.00000000e+00,
1.83424060e-20, 0.00000000e+00, -1.83424060e-20, 0.00000000e+00,
9.17120299e-21, 0.00000000e+00, -1.83424060e-21]),
array([ 1. , -9.99851389, 44.98765092, -119.95470631,
209.90388506, -251.87018009, 209.88453023, -119.93258575,
44.9752074 , -9.99482662, 0.99953904]))
查看b
系数(第一个数组):它们的值为 1e-20,这意味着滤波器设计完全无法收敛,如果将其应用于任何信号,输出将为零——这就是你发现的。
您没有提及您的应用程序,但如果您真的想将信号的频率内容保持在 100 到 101 Hz 之间,您可以对信号进行零填充 FFT,将该频带之外的频谱部分归零,然后 IFFT (查看rfft
、irfft
和rfftfreq
在numpy.fft
模块中)。
这是一个使用 FFT 在傅里叶域中应用砖墙带通滤波器的函数:
import numpy.fft as fft
import numpy as np
def fftBandpass(x, low, high, fs=1.0):
"""
Apply a bandpass signal via FFTs.
Parameters
----------
x : array_like
Input signal vector. Assumed to be real-only.
low : float
Lower bound of the passband in Hertz. (If less than or equal
to zero, a high-pass filter is applied.)
high : float
Upper bound of the passband, Hertz.
fs : float
Sample rate in units of samples per second. If `high > fs / 2`,
the output is low-pass filtered.
Returns
-------
y : ndarray
Output signal vector with all frequencies outside the `[low, high]`
passband zeroed.
Caveat
------
Note that the energe in `y` will be lower than the energy in `x`, i.e.,
`sum(abs(y)) < sum(abs(x))`.
"""
xf = fft.rfft(x)
f = fft.rfftfreq(len(x), d=1 / fs)
xf[f < low] = 0
xf[f > high] = 0
return fft.irfft(xf, len(x))
if __name__ == '__main__':
fs = 44.1e3
N = int(fs)
x = np.random.randn(N)
t = np.arange(N) / fs
import pylab as plt
plt.figure()
plt.plot(t, x, t, 100 * fftBandpass(x, 100, 101, fs=fs))
plt.xlabel('time (seconds)')
plt.ylabel('signal')
plt.legend(['original', 'scaled bandpassed'])
plt.show()
您可以将它放在一个文件中fftBandpass.py
,然后运行它python fftBandpass.py
以查看它创建以下图:
请注意,我必须将 1 Hz 带通信号缩放 100,因为在带通那么多之后,信号中的能量非常少。另请注意,位于这个小通带内的信号几乎只是 100 Hz 左右的正弦曲线。
如果您将以下内容放入您自己的代码中:from fftBandpass import fftBandpass
,您可以使用该fftBandpass
功能。
您可以尝试的另一件事是将信号抽取 100 倍,因此将其转换为以 441 Hz 采样的信号。441 Hz 中的 1 Hz 仍然是一个非常窄的通带,但您可能比尝试带通原始信号有更好的运气。请参阅scipy.signal.decimate
,但不要尝试使用 调用它q=100
,而是递归抽取信号,依次为 2、2、5、5(总抽取 100 倍)。
推荐阅读
- google-sheets - 根据某些特点选择最佳方案
- angular - JSON.stringify之后的Angular4材料日期选择器时区问题
- ios - 为应用程序 ID 启用 iCloud Storage API
- jquery - m.自定义滚动条在 RTL 模式下不起作用
- c# - 初始化匿名类型实体框架时如何使用虚拟属性
- java - Apache Cayenne - 批量关系设置
- r - 无法使用 remove() 删除 r 中的变量
- lua - ZeroBrane - 断点未命中
- google-cloud-platform - 在 Cloud SQL 中批量加载数据的 ETL 方法
- f# - 使用 FParsec 解析多个语句时不显示错误