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

Question

如何应用低通滤波器，截止频率线性变化（或曲线比线性更一般），例如10000hz 到 200hz 的时间，使用 numpy/scipy，可能没有其他库？

例子：

• 在 00:00,000，低通截止 = 10000hz
• 在 00:05,000，低通截止 = 5000hz
• 在 00:09,000，低通截止 = 1000hz
• 然后截止频率在 10 秒内保持在 1000hz，然后截止频率降低到 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))

但是如何使截止值变化？

注意：我已经阅读了Applying time-variant filter in Python，但答案相当复杂（通常适用于多种过滤器）。

Answer 1

一种相对简单的方法是保持滤波器固定并改为调制信号时间。例如，如果信号时间快 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 秒的频谱图，通过时间弯曲实现了时变低通滤波。

Answer 2

您可以使用 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 赫兹的频率设置为零

对于时变截止，我会拆分声音并应用过滤器。假设声音的采样率为 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）