python - 在Python中应用时变过滤器

我正在尝试使用一个带时变截止频率的带通滤波器来使用Python。我目前使用的例程将我的信号分割成等长度的时间段，然后在每次将信号合并之前，对每个段应用一个具有特定于时间的参数的滤波器。参数是基于预先存在的估计。我所遇到的问题是，在应用了过滤器之后出现的每个时间段的边缘都有“涟漪”。这会导致我的信号不连续，干扰我的滤波后数据分析。我希望有人能告诉我，在Python中是否存在用于实现具有时变参数的过滤器的现有例程？另外，关于我如何解决这个问题的建议将非常感谢。编辑——下面添加了我想做的事情的示例。假设我有一个信号x（t）。我想用参数（100200）Hz的带通滤波器对上半部分进行滤波。下半部分我想用参数（140, 240）Hz滤波。我迭代X（t），将我的过滤器应用到每一半，然后重新组合结果。一些示例代码可能如下所示：outputArray = np.empty(len(x))segmentSize = len(x) / 2filtParams = [(100, 200), (140, 240)]for i in range(2): tempData = x[i*segmentSize:(i+1)*segmentSize] tempFiltered = bandPassFilter(tempData, parameters=filtParams[i]) outputArray[i*segmentSize:(i+1)*segmentSize] = tempFiltered（为了节省空间，假设我有一个执行带通滤波的函数）。如您所见，数据段不重叠，只是简单地“粘贴”在新数组中。编辑2——我的问题@ HD的一些示例代码。首先，感谢你迄今为止的重要投入。AudioLasy软件包看起来是个很棒的工具。我想如果我能更详细地描述我的目标会更有用一些。正如我已经发表的elsewhere，我试图用希尔伯特变换提取一个信号的AA>（if）。我的数据包含了很大的噪声，但是我对我的IF信号的带宽有很好的估计。不过，我遇到的一个问题是IF常常是非平稳的。因此，使用“静态”滤波器的方法经常需要使用宽的带通区域，以确保捕获所有频率。下面的代码演示了增加滤波器带宽对中频信号的影响。它包括信号生成功能、使用SIPY.Stand信号包的带通滤波器的实现以及提取所得到的滤波信号的IF的方法。from audiolazy import *import scipy.signal as sigimport numpy as npfrom pylab import *def sineGenerator( ts, f, rate, noiseLevel=None ): """generate a sine tone with time, frequency, sample rate and noise specified""" fs = np.ones(len(ts)) * f y = np.sin(2*np.pi*fs*ts) if noiseLevel: y = y + np.random.randn(len(y))/float(noiseLevel) return ydef bandPassFilter( y, passFreqs, rate, order ): """STATIC bandpass filter using scipy.signal Butterworth filter""" nyquist = rate / 2.0 Wn = np.array([passFreqs[0]/nyquist, passFreqs[1]/nyquist]) z, p, k = sig.butter(order, Wn, btype='bandpass', output='zpk') b, a = sig.zpk2tf(z, p, k) return sig.lfilter(b, a, y)if __name__ == '__main__': rate = 1e4 ts = np.arange(0, 10, 1/rate) # CHANGING THE FILTER AFFECTS THE LEVEL OF NOISE ys = sineGenerator(ts, 600.0, 1e4, noiseLevel=1.0) # a 600Hz signal with noise filts = [[500, 700], [550, 650], [580, 620]] for f in filts: tempFilt = bandPassFilter( ys, f, rate, order=2 ) tempFreq = instantaneousFrequency( tempFilt, rate ) plot( ts[1:], tempFreq, alpha=.7, label=str(f).strip('[]') ) ylim( 500, 750 ) xlabel( 'time' ) ylabel( 'instantaneous frequency (Hz)' ) legend(frameon=False) title('changing filter passband and instantaneous frequency') savefig('changingPassBand.png')信号中只有一个频率分量（600Hz）。图例显示了在每种情况下使用的滤波器参数。使用一个更窄的“静态”过滤器给出一个“干净”的输出。但是我的滤波器能有多窄受到频率的限制。例如，考虑以下两个频率分量的信号（一个在600Hz，另一个在650Hz）。在这个例子中，我被迫使用一个更宽的带通滤波器，这导致额外的噪声悄悄进入中频数据。我的想法是，通过使用时变滤波器，我可以在某些时间增量上“优化”我的信号的滤波器。所以对于上面的信号，我可能希望在前5秒过滤580-620Hz，然后在后5秒过滤630-670Hz。基本上我想在最终的中频信号中最小化噪声。基于您发送的示例代码，我编写了一个函数，该函数使用AdoLoZy在信号上实现静态巴特沃思滤波器。def audioLazyFilter( y, rate, Wp, Ws ): """implement a Butterworth filter using audiolazy""" s, Hz = sHz(rate) Wp = Wp * Hz # Bandpass range in rad/sample Ws = Ws * Hz # Bandstop range in rad/sample order, new_wp_divpi = sig.buttord(Wp/np.pi, Ws/np.pi, gpass=dB10(.6), gstop=dB10(.1)) ssfilt = sig.butter(order, new_wp_divpi, btype='bandpass') filt_butter = ZFilter(ssfilt[0].tolist(), ssfilt[1].tolist()) return list(filt_butter(y))使用该滤波器结合希尔伯特变换例程获得的IF数据与使用SiPi.Stand获得的那些信号进行了比较：AL_filtered = audioLazyFilter( ys, rate, np.array([580, 620]), np.array([570, 630]) )SP_filtered = bandPassFilter( ys, [580, 620], rate, order=2 )plot(ts[1:], instantaneousFrequency( SP_filtered, 1e4 ), alpha=.75, label='scipy.signal Butterworth filter')plot(ts[1:], instantaneousFrequency( AL_filtered, 1e4 ), 'r', alpha=.75, label='audiolazy Butterworth filter')我的问题是现在我可以将AdioLasZy ButtValm例程与您在原始文章中描述的时变特性结合起来吗？ 最佳答案 使用时变滤波器的本地工作原理from audiolazy import sHz, white_noise, line, resonator, AudioIOrate = 44100s, Hz = sHz(rate)sig = white_noise() # Endless white noise Streamdur = 8 * s # Some few seconds of audiofreq = line(dur, 200, 800) # A lazy iterable rangebw = line(dur, 100, 240)filt = resonator(freq * Hz, bw * Hz) # A simple bandpass filterwith AudioIO(True) as player: player.play(filt(sig), rate=rate)您也可以使用它来绘制（或一般的）分析，通过使用list(filt(sig))或 >。还有很多其他的资源可能是有用的，例如在z变换滤波器方程中直接应用时变系数。编辑：有关AudioLazy组件的详细信息下面的例子是使用ipython完成的。谐振器是一个cc实例，它将许多实现绑定在一个地方。In [1]: from audiolazy import *In [2]: resonatorOut[2]:{('freq_poles_exp',): <function audiolazy.lazy_filters.freq_poles_exp>, ('freq_z_exp',): <function audiolazy.lazy_filters.freq_z_exp>, ('poles_exp',): <function audiolazy.lazy_filters.poles_exp>, ('z_exp',): <function audiolazy.lazy_filters.z_exp>}In [3]: resonator.defaultOut[3]: <function audiolazy.lazy_filters.poles_exp>因此filt(sig).take(inf)在内部调用StrategyDict函数，从中可以获得一些帮助In [4]: resonator.poles_exp?Type: functionString Form:<function poles_exp at 0x2a55b18>File: /usr/lib/python2.7/site-packages/audiolazy/lazy_filters.pyDefinition: resonator.poles_exp(freq, bandwidth)Docstring:Resonator filter with 2-poles (conjugated pair) and no zeros (constantnumerator), with exponential approximation for bandwidth calculation.Parameters----------freq : Resonant frequency in rad/sample (max gain).bandwidth : Bandwidth frequency range in rad/sample following the equation: ``R = exp(-bandwidth / 2)`` where R is the pole amplitude (radius).Returns-------A ZFilter object.Gain is normalized to have peak with 0 dB (1.0 amplitude).所以一个详细的过滤器分配将是filt = resonator.poles_exp(freq=freq * Hz, bandwidth=bw * Hz)其中resonator只是一个数字，将单元从HZ转换到RAAD/SAMPLE，就像大多数音频文件中使用的一样。让我们使用resonator.poles_exp和Hz（freq = pi/4已经在bw = pi/8命名空间中）：In [5]: filt = resonator(freq=pi/4, bandwidth=pi/8)In [6]: filtOut[6]: 0.233921------------------------------------1 - 1.14005 * z^-1 + 0.675232 * z^-2In [7]: type(filt)Out[7]: audiolazy.lazy_filters.ZFilter您可以尝试使用此筛选器，而不是第一个示例中给定的筛选器。另一种方法是使用包中的pi对象。首先让我们找出全极点谐振器的常数：In [8]: freq, bw = pi/4, pi/8In [9]: R = e ** (-bw / 2)In [10]: c = cos(freq) * 2 * R / (1 + R ** 2) # AudioLazy included the cosineIn [11]: gain = (1 - R ** 2) * sqrt(1 - c ** 2)分母可以直接使用公式中的audiolazy来完成：In [12]: denominator = 1 - 2 * R * c * z ** -1 + R ** 2 * z ** -2In [13]: gain / denominatorOut[14]: 0.233921------------------------------------1 - 1.14005 * z^-1 + 0.675232 * z^-2In [15]: type(_) # The "_" is the last returned value in IPythonOut[15]: audiolazy.lazy_filters.ZFilter编辑2：关于时变系数滤波器系数也可以是流实例（可以从任何iterable转换）。In [16]: coeff = Stream([1, -1, 1, -1, 1, -1, 1, -1, 1, -1]) # Cast from a listIn [17]: (1 - coeff * z ** -2)(impulse()).take(inf)Out[17]: [1.0, 0.0, -1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]同样，给出一个列表输入而不是z流：In [18]: coeff = Stream((1, -1, 1, -1, 1, -1, 1, -1, 1, -1)) # Cast from a tupleIn [19]: (1 - coeff * z ** -2)([1, 0, 0, 0, 0, 0, 0]).take(inf)Out[19]: [1.0, 0.0, -1, 0, 0, 0, 0]numpy 1d数组也是一个iterable：In [20]: from numpy import arrayIn [21]: array_data = array([1, -1, 1, -1, 1, -1, 1, -1, 1, -1])In [22]: coeff = Stream(array_data) # Cast from an arrayIn [23]: (1 - coeff * z ** -2)([0, 1, 0, 0, 0, 0, 0]).take(inf)Out[23]: [0.0, 1.0, 0, 1, 0, 0, 0]最后一个例子展示了时变行为。编辑3：分块重复序列行为line函数的行为类似于z，它获取范围“length”而不是“step”。In [24]: import numpyIn [25]: numpy.linspace(10, 20, 5) # Start, stop (included), lengthOut[25]: array([ 10. , 12.5, 15. , 17.5, 20. ])In [26]: numpy.linspace(10, 20, 5, endpoint=False) # Makes stop not includedOut[26]: array([ 10., 12., 14., 16., 18.])In [27]: line(5, 10, 20).take(inf) # Length, start, stop (range-like)Out[27]: [10.0, 12.0, 14.0, 16.0, 18.0]In [28]: line(5, 10, 20, finish=True).take(inf) # Include the "stop"Out[28]: [10.0, 12.5, 15.0, 17.5, 20.0]这样，过滤器方程每个样本具有不同的行为样本（1样本“块”）。不管怎样，对于较大的块大小，可以使用中继器：In [29]: five_items = _ # List from the last Out[] valueIn [30]: @tostream ....: def repeater(sig, n): ....: for el in sig: ....: for _ in xrange(n): ....: yield el ....:In [31]: repeater(five_items, 2).take(inf)Out[31]: [10.0, 10.0, 12.5, 12.5, 15.0, 15.0, 17.5, 17.5, 20.0, 20.0]并在第一个示例的行中使用它，使impulse()和numpy.linspace变成：chunk_size = 100freq = repeater(line(dur / chunk_size, 200, 800), chunk_size)bw = repeater(line(dur / chunk_size, 100, 240), chunk_size)编辑4：利用时变增益/包络仿真LTI滤波器的时变滤波器/系数另一种方法是使用不同的“权重”来对两个不同的滤波版本进行信号处理，并用信号进行一些“交叉渐变”的数学运算，比如：signal = thub(sig, 2) # T-Hub is a T (tee) auto-copyfilt1(signal) * line(dur, 0, 1) + filt2(signal) * line(dur, 1, 0)这将从同一信号的不同滤波版本应用线性包络（从0到1和从1到0）。如果 >看起来很混乱，尝试使用“cc”和“cc”，而这些应该做同样的事情。编辑5：时变Butterworth滤波器我花了最后几个小时试图让北海个人化为你的例子，强加freq并给予半功率带宽（~3dB）直接。我已经做了四个例子，代码是AA>，并且我已经更新了AdioLoZy以包含一个bw高斯分布的噪声流。请注意，GIST中的代码没有优化，它是在这个特定的情况下工作的，并且Chirp示例使得它由于“每个样本”系数查找行为而变得非常慢。瞬时频率可以从/采样中的[过滤]数据得到：diff(unwrap(phase(hilbert(filtered_data))))在thub或另一种在离散时间内找到导数的方法中，sig1, sig2 = tee(sig, 2)是从“cc>函数中给出的解析信号（离散希尔伯特变换是它的结果的虚部），另外两个是来自AdioLoZy的辅助函数。这是当巴特沃思突然改变其系数，同时保持其记忆，没有噪声：在这种转变中，一种明显的振荡行为是显而易见的。你可以使用一个移动中值来“平滑”在低频侧，同时保持突然的转变，但是这不会用更高的频率来工作。这是我们期望的完美正弦曲线，但是有噪声（很多噪声，高斯的标准偏差等于正弦振幅），它变成：我试着用同样的声音来做，正是这样：这显示了一个奇怪的行为时，低带宽滤波，在最高频率。加性噪声：gist中的代码也filt(sig1)这最后一个噪声啁啾。编辑6：时变谐振器滤波器我在AudioLazy中添加了一个使用谐振器而不是巴特沃斯的例子。它们在纯Python中并没有被优化，但是在啁啾中对每个样本的调用速度比调用filt(sig2)要快，并且更容易实现，因为所有order = 2策略都接受流实例作为有效输入。以下是两个谐振器级联（即二阶滤波器）的曲线图：对于三个谐振器的级联（即，第三阶滤波器）是相同的：这些谐振器在中心频率处的增益等于1（0db），并且即使没有任何滤波，从过渡过程中的“突变纯正弦波”图中的振荡模式也会发生。关于python - 在Python中应用时变过滤器，我们在Stack Overflow上找到一个类似的问题：https://stackoverflow.com/questions/18128057/