具有时间变化截止频率的低通滤波器，使用Python实现

Question

具有时间变化截止频率的低通滤波器，使用Python实现

11

如何使用numpy/scipy应用低通滤波器，并使截止频率沿时间线性变化（或者采用比线性更一般的曲线），例如从10000Hz到200Hz？

示例：

- 在00:00,000时，低通截止频率为10000Hz - 在00:05,000时，低通截止频率为5000Hz - 在00:09,000时，低通截止频率为1000Hz - 然后截止频率保持在1000Hz，持续10秒钟，然后下降到200Hz

以下是如何实现简单的100Hz低通滤波器的方法：

from scipy.io import wavfile
import numpy as np
from scipy.signal import butter, lfilter

sr, x = wavfile.read('test.wav')
b, a = butter(2, 100.0 / sr, btype='low')  # Butterworth
y = lfilter(b, a, x)
wavfile.write('out.wav', sr, np.asarray(y, dtype=np.int16))

但如何使截止频率变化？

注意：我已经阅读了在Python中应用时变滤波器，但答案相当复杂（而且通常适用于许多类型的滤波器）。

- Basj

总是二阶巴特沃斯滤波器吗？ - Stephen Rauch

@StephenRauch 任何IIR或FIR都可以，只要截止频率C（t）可以沿着时间t平稳地变化（即C（t）不是分段常数函数）。 - Basj

听起来像是需要使用小波变换，但我对它们的经验不足以给出一个好的答案。不幸的是，pywavelets和pywt似乎已经过时了。 - Daniel F

初始采样率是多少？ - Daniel F

@DanielF 44.1千赫、48千赫或96千赫是这种情况下常用的典型值。 - Basj

2个回答

3

您可以使用scipy.fftpack.fftfreq和scipy.fftpack.rfft来设置阈值。

fft = scipy.fftpack.fft(sound)
freqs = scipy.fftpack.fftfreq(sound.size, time_step)

对于time_step，我使用了音频采样率的两倍

fft[(freqs < 200)] = 0

这将把小于200 hz的所有频率都设为零。

对于随时间变化的截止频率，我会先拆分音频，再应用滤波器。假设音频采样率为44100，则5000Hz的滤波器将从第220500个样本（即5秒）开始。

10ksound = sound[:220500]
10kfreq = scipy.fftpack.fftreq(10ksound.size, time_step)
10kfft = scipy.fftpack.fft(10ksound)
10kfft[(10kfreqs < 10000)] = 0

接下来是下一个过滤器：

5ksound = sound[220500:396900]
5kfreq = scipy.fftpack.fftreq(10ksound.size, time_step)
5kfft = scipy.fftpack.fft(10ksound)
5kfft[(5kfreqs < 5000)] = 0

编辑：为了使其“滑动”或者渐变，你可以将“块”变得更小，并且对应的每个块应用越来越大的频率阈值(5000 -> 5001 -> 5002)。

- Jacob Sharp

谢谢你的答复，有几点需要说明：1）将 bins 归零不太适合用于过滤（但它有点有效），请参阅有关 DSP 的主题，2）低通意味着我们保留低频并消除高频（您已经相反了，这很容易修复）。最大的问题是：3）实际上您的过滤器是声音每个部分的 2 或 3 个标准过滤器。您具有“分段常数截止”，这很简单。难点在于具有“滑动”截止（请参见问题）。 - Basj

你是正确的，对于小块（例如1024个样本，即在44.1Khz采样率下约为23ms），可以在“滑动情况”下工作：https://dev59.com/aa7la4cB1Zd3GeqPfpmK#52469397。重要的部分是关注初始条件参数。 - Basj

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- Paul Panzer · Accepted Answer

一种相对简单的方法是保持过滤器固定，并改变信号时间。例如，如果信号时间运行速度比标准时间快10倍，则10kHz低通滤波器将起到1kHz低通滤波器的作用。

为了实现这一点，我们需要解决一个简单的ODE。

dy       1
--  =  ----
dt     f(y)

这里的t是调制时间，y是实际时间，f是在时间y时所需的截止频率。

原型实现：

from __future__ import division
import numpy as np
from scipy import integrate, interpolate
from scipy.signal import butter, lfilter, spectrogram

slack_l, slack = 0.1, 1
cutoff = 50
L = 25

from scipy.io import wavfile
sr, x = wavfile.read('capriccio.wav')
x = x[:(L + slack) * sr, 0]
x = x

# sr = 44100
# x = np.random.normal(size=((L + slack) * sr,))

b, a = butter(2, 2 * cutoff / sr, btype='low')  # Butterworth

# cutoff function
def f(t):
    return (10000 - 1000 * np.clip(t, 0, 9) - 1000 * np.clip(t-19, 0, 0.8)) \
        / cutoff

# and its reciprocal
def fr(_, t):
    return cutoff / (10000 - 1000 * t.clip(0, 9) - 1000 * (t-19).clip(0, 0.8))

# modulate time
# calculate upper end of td first
tdmax = integrate.quad(f, 0, L + slack_l, points=[9, 19, 19.8])[0]
span = (0, tdmax)
t = np.arange(x.size) / sr
tdinfo = integrate.solve_ivp(fr, span, np.zeros((1,)),
                             t_eval=np.arange(0, span[-1], 1 / sr),
                             vectorized=True)
td = tdinfo.y.ravel()
# modulate signal
xd = interpolate.interp1d(t, x)(td)
# and linearly filter
yd = lfilter(b, a, xd)
# modulate signal back to linear time
y = interpolate.interp1d(td, yd)(t[:-sr*slack])

# check
import pylab
xa, ya, z = spectrogram(y, sr)
pylab.pcolor(ya, xa, z, vmax=2**8, cmap='nipy_spectral')
pylab.savefig('tst.png')

wavfile.write('capriccio_vandalized.wav', sr, y.astype(np.int16))

示例输出：

BWV 826 Capriccio前25秒的频谱图被通过时间弯曲实现的时变低通滤波器进行了过滤。