在Python中应用时变滤波器

Question 1

我试图用Python对一个信号应用一个具有时间变化的截止频率的带通滤波器。我目前使用的程序是将我的信号分成等长的时间段，然后对每个时间段应用一个具有特定时间参数的滤波器，然后再将信号合并起来。这些参数是基于预先存在的估计。

我遇到的问题是，在每个时间段的边缘有 "波纹"，是在应用滤波器后出现的。这导致我的信号出现不连续，干扰了我的滤波后数据分析。

我希望有人能告诉我，在Python中是否有任何现有的程序来实现具有时变参数的过滤器？另外，如果有关于如何解决这个问题的建议，我将非常感激。

EDIT -- 我想做的事情的例子补充如下。

假设我有一个信号x(t)。我想用一个参数为(100,200)Hz的带通滤波器来过滤前半部分。后半部分我想用参数（140，240）Hz来过滤。我对x(t)进行迭代，将我的滤波器应用于每一半，然后将结果重新组合。一些示例代码可能看起来像。

outputArray = np.empty(len(x))
segmentSize = len(x) / 2
filtParams  = [(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
(为了节省篇幅，我们假设我有一个进行带通滤波的函数）。
正如你所看到的，这些数据段没有重叠，只是被 "粘贴 "回新的数组中。
EDIT 2-- 我的问题的一些示例代码 @H.D.
首先，感谢你到目前为止的重要投入。audiolazy软件包看起来是一个伟大的工具。
我想如果我进一步详细描述我的目标，会更有用一些。正如我已经发布的elsewhere，我正试图提取瞬时频率(使用希尔伯特变换，对一个信号的中频（IF）进行测量。我的数据含有大量的噪声，但我对中频信号所在的带宽有一个很好的估计。然而，我遇到的一个问题是
中频往往是非稳态的。因此，使用 "静态 "滤波器的方法，我经常需要使用一个宽的带通区域，以确保所有的频率都被捕获。
下面的代码演示了增加滤波器带宽对中频信号的影响。它包括一个信号生成函数，一个使用scipy.signal包的带通滤波器的实现，以及一个提取所得滤波信号中频的方法。
from audiolazy import *
import scipy.signal as sig
import numpy as np
from 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 y
def 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')
在这个例子中，我被迫使用了一个更宽的带通滤波器，这导致了额外的噪声悄悄进入中频数据。
我的想法是，通过使用一个时间变化的滤波器，我可以在某些时间增量上为我的信号 "优化 "滤波器。因此，对于上述信号，我可能想在前5秒过滤580-620Hz，然后在后5秒过滤630-670Hz。从本质上讲，我想尽量减少最终中频信号的噪音。
根据你发送的示例代码，我写了一个函数，使用audiolazy在信号上实现一个静态巴特沃斯滤波器。
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))
使用这个滤波器与希尔伯特变换例程一起获得的中频数据与使用scipy.signal获得的数据相比，效果很好。
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')
我的问题是，现在我可以将audiolazy Butterworth程序与你在原帖中描述的时变特性结合起来吗？

Question 2


          
           
            
             
              
               懒惰的音频
              
              与时间变化的过滤器一起工作
             
             from audiolazy import sHz, white_noise, line, resonator, AudioIO
rate = 44100
s, Hz = sHz(rate)
sig = white_noise() # Endless white noise Stream
dur = 8 * s # Some few seconds of audio
freq = line(dur, 200, 800) # A lazy iterable range
bw = line(dur, 100, 240)
filt = resonator(freq * Hz, bw * Hz) # A simple bandpass filter
with AudioIO(True) as player:
  player.play(filt(sig), rate=rate)
你也可以通过使用list(filt(sig))或filt(sig).take(inf)，将其用于绘图（或分析，一般来说）。还有很多其他的资源可能也是有用的，比如在Z-变换滤波器方程中直接应用时变系数。
EDIT: More information about the 懒惰的音频 components
下面的例子是用IPython完成的。
共振器是一个StrategyDict的实例，它将许多实现绑在一个地方。
In [1]: from audiolazy import *
In [2]: resonator
Out[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.default
Out[3]: <function audiolazy.lazy_filters.poles_exp>
所以resonator在内部调用resonator.poles_exp函数，你可以从中得到一些帮助
In [4]: resonator.poles_exp?
Type:       function
String Form:<function poles_exp at 0x2a55b18>
File:       /usr/lib/python2.7/site-packages/audiolazy/lazy_filters.py
Definition: resonator.poles_exp(freq, bandwidth)
Docstring:
Resonator filter with 2-poles (conjugated pair) and no zeros (constant
numerator), 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)
Where the Hz is just a number to change the unit from Hz to rad/sample, as used in most 懒惰的音频 components.
Let's do so with freq = pi/4 and bw = pi/8 (pi is already in the audiolazy namespace):
In [5]: filt = resonator(freq=pi/4, bandwidth=pi/8)
In [6]: filt
Out[6]: 
              0.233921
------------------------------------
1 - 1.14005 * z^-1 + 0.675232 * z^-2
In [7]: type(filt)
Out[7]: audiolazy.lazy_filters.ZFilter
你可以尝试使用这个过滤器，而不是第一个例子中给出的那个。
另一种方法是使用软件包中的z对象。
首先让我们找到那个全极谐振器的常数。
In [8]: freq, bw = pi/4, pi/8
In [9]: R = e ** (-bw / 2)
In [10]: c = cos(freq) * 2 * R / (1 + R ** 2) # AudioLazy included the cosine
In [11]: gain = (1 - R ** 2) * sqrt(1 - c ** 2)
分母可以直接用方程中的z来完成。
In [12]: denominator = 1 - 2 * R * c * z ** -1 + R ** 2 * z ** -2
In [13]: gain / denominator
Out[14]: 
              0.233921
------------------------------------
1 - 1.14005 * z^-1 + 0.675232 * z^-2
In [15]: type(_) # The "_" is the last returned value in IPython
Out[15]: audiolazy.lazy_filters.ZFilter
编辑2：关于时间变化的系数
滤波器的系数也可以是一个Stream实例（可以从任何迭代器中投出）。
In [16]: coeff = Stream([1, -1, 1, -1, 1, -1, 1, -1, 1, -1]) # Cast from a list
In [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]
同样，给定一个列表输入而不是impulse()。流。
In [18]: coeff = Stream((1, -1, 1, -1, 1, -1, 1, -1, 1, -1)) # Cast from a tuple 
In [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数组也是一个可迭代的。
In [20]: from numpy import array
In [21]: array_data = array([1, -1, 1, -1, 1, -1, 1, -1, 1, -1])
In [22]: coeff = Stream(array_data) # Cast from an array
In [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]
这最后一个例子显示了时间变化的行为。
EDIT 3:分块重复的序列行为
线条函数的行为类似于numpy.linspace，它得到的范围是 "长度 "而不是 "步长"。
In [24]: import numpy
In [25]: numpy.linspace(10, 20, 5) # Start, stop (included), length
Out[25]: array([ 10. ,  12.5,  15. ,  17.5,  20. ])
In [26]: numpy.linspace(10, 20, 5, endpoint=False) # Makes stop not included
Out[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[] value
In [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]
并在第一个例子的行中使用，使freq和bw成为。
chunk_size = 100
freq = 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-copy
filt1(signal) * line(dur, 0, 1) + filt2(signal) * line(dur, 1, 0)
这将从同一信号的不同过滤版本中应用一个线性包络（从0到1和从1到0）。如果thub看起来很混乱，请尝试sig1, sig2 = tee(sig, 2)应用filt(sig1)和filt(sig2)，这些应该做同样的事情。
EDIT 5：时变巴特沃斯滤波器
我花了几个小时试图让那个巴特沃斯像你的例子一样个性化，施加order = 2并直接给出半功率带宽（~3dB）。我已经做了四个例子，代码是in this Gist, and I've updated 懒惰的音频 to include a gauss_noise Gaussian-distributed noise stream. Please note that the code in gist has nothing optimized, it was done ony to work in this particular case, and the chirp example makes it really slow due to a "per sample" coefficient finding behaviour. The instant frequency can be get from the [filtered] data in rad/sample with:
diff(unwrap(phase(hilbert(filtered_data))))
where diff = 1 - z ** -1 or another approach to find derivatives in discrete time, hilbert is the function from scipy.signal that gives us the analytical signal (the Discrete Hilbert Transform is the imaginary part of its result) and the other two are helper functions from 懒惰的音频.
这就是巴特沃斯在保持记忆的同时突然改变其系数的情况，没有噪音。
在这个过渡中，明显有振荡行为。你可以用移动中值来 "平滑 "低频部分，同时保持突然的过渡，但这对高频部分不起作用。好吧，这就是我们从一个完美的正弦波中所期望的，但如果有了噪声（大量的噪声，高斯的标准差等于正弦波的振幅），它就变成了。
我当时试图用鸣笛做同样的事情，正是这样。
当用较低的带宽进行滤波时，在最高频率上显示出一种奇怪的行为。而且是在有加法噪声的情况下。
gist中的代码也AudioIO().play这个最后的嘈杂的鸣叫。
EDIT 6: 时变谐振器滤波器
我已经添加到大意如此一个使用谐振器而不是巴特沃斯的例子。它们是纯Python语言，没有经过优化，但比在啁啾过程中对每个样本调用butter要快，而且更容易实现，因为所有的resonator策略都接受Stream实例作为有效输入。下面是两个谐振器级联的图示（即二阶滤波器）。
这些谐振器在中心频率上的增益等于1（0dB），即使没有任何滤波，也会发生 "突然的纯正弦波 "过渡图中的那种振荡模式。

Question 3


          
           
            
             
              
               
                
                 
                  
                   
                    
                     
                      
                       
                        
                         
                          
                           
                            
                             试着用你所使用的每个滤波器过滤整个信号，然后适当地合并过滤后的信号。粗略的说，用伪代码。
                            
                            # each filter acts on the entire signal with each cutoff frequency
s1 = filter1(signal)
s2 = filter2(signal)
s3 = filter3(signal)
filtered_signal = s1[0:t1] + s2[t1:t2] + s3[t2:t3]
我认为这将避免你所描述的因切碎信号然后进行过滤而产生的伪影。
另一种可能性是使用一个短时傅里叶变换(STFT).  Here's 使用numpy的实现.
基本上，你可以对你的信号进行STFT，通过对时间-频率阵列的操作来过滤你的信号，然后反STFT阵列来得到一个过滤的信号。
你可能会得到更好的结果，使用一个可倒置的小波变换.  这里有一个屏蔽的论文 describing how to do pretty much what you're trying to do with a 小波变换.

Question 4


          
           
            
             
              
               
                
                 
                  
                   
                    
                     
                      
                       
                        
                         
                          
                           
                            
                             
                              如果你用切片提取信号的一部分，那么你实际上是用一个矩形窗口来给你的数据开窗，由于突然的不连续，在边缘会有 "响声"。解决这个问题的一个方法是使用一个响声较小的窗口，如汉宁窗口。
                             
                             import numpy as np
signal = np.random.randn(222)
length = 50
window = np.hanning(length)