具有时变截止频率的低通滤波器，使用Python [英] Lowpass filter with a time-varying cutoff frequency, with Python

问题描述

如何应用低通滤波器，其截止频率从例如随时间变化从10000hz到200hz ，有numpy/scipy，可能没有其他库了?

How to apply a lowpass filter, with cutoff frequency varying linearly (or with a more general curve than linear) from e.g. 10000hz to 200hz along time, with numpy/scipy and possibly no other library?

示例:

• 在00:00,000时，低通截止频率= 10000hz
• 在00:05,000时，低通截止频率= 5000hz
• 00:09,000，低通截止频率= 1000hz
• 然后截止频率在10秒内保持在1000hz，然后截止频率降低到200hz

这里是执行简单的100hz低通的方法:

Here is how to do a simple 100hz lowpass:

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))

但是如何使临界值变化?

but how to make the cutoff vary?

注意:我已经阅读了在Python中应用时变过滤器，但答案非常复杂(通常适用于多种过滤器).

Note: I've already read Applying time-variant filter in Python but the answer is quite complex (and it applies to many kinds of filter in general).

推荐答案

一种比较简单的方法是保持滤波器固定不变，而是调制信号时间.例如，如果信号时间快10倍，则在标准时间内10KHz低通将像1KHz低通一样.

One comparatively easy method is to keep the filter fixed and modulate signal time instead. For example, if signal time runs 10x faster a 10KHz lowpass will act like a 1KHz lowpass in standard time.

为此，我们需要解决一个简单的ODE

To do this we need to solve a simple ODE

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

此处t是调制时间y的实时时间，f是时间y的所需截止时间.

Here t is modulated time y real time and f the desired cutoff at time 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秒频谱图，通过通过时间弯曲实现的时变低通进行了滤波.

Spectrogram of first 25 seconds of BWV 826 Capriccio filtered with a time varying lowpass implemented via time bending.

这篇关于具有时变截止频率的低通滤波器，使用Python的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！