5

我正在尝试使用 Python 将具有时变截止频率的带通滤波器应用于信号。我目前使用的例程将我的信号划分为等长的时间段,然后对于每个段,我应用具有特定时间参数的过滤器,然后将信号重新合并在一起。这些参数基于预先存在的估计。

我似乎遇到的问题是在应用过滤器后出现的每个时间段的边缘都有“涟漪”。这会导致我的信号不连续,从而干扰我的滤波后数据分析。

我希望有人能告诉我是否有任何现有的例程可以在 Python 中实现具有时变参数的过滤器?或者,我将不胜感激有关如何解决此问题的建议。

编辑——下面添加了我想要做的示例。

假设我有一个信号 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

(为了节省空间,假设我有一个执行带通滤波的函数)。

如您所见,数据段不重叠,而是简单地“粘贴”回新数组中。

编辑 2——我的问题的一些示例代码@HD

首先,感谢您迄今为止的重要投入。audiolazy 包看起来是个很棒的工具。

我认为如果我更详细地描述我的目标会更有用。正如我在其他地方发布的那样,我正在尝试使用希尔伯特变换来提取信号的瞬时频率(IF)。我的数据包含大量噪声,但我对 IF 信号所在的带宽有一个很好的估计。然而,我遇到的一个问题是 IF 通常是非平稳的。因此,使用“静态”滤波器方法,我经常需要使用宽带通区域,以确保捕获所有频率。

以下代码演示了增加滤波器带宽对 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')

改变通带

信号中有一个单一的频率分量(600Hz)。图例显示了每种情况下使用的过滤器参数。使用更窄的“静态”过滤器可以提供“更干净”的输出。但是我的滤波器可以有多窄受频率的限制。例如,考虑以下具有两个频率分量的信号(一个在 600Hz,另一个在 650Hz)。

变频

在此示例中,我被迫使用更宽的带通滤波器,这导致额外的噪声潜入 IF 数据。

我的想法是,通过使用时变滤波器,我可以在特定时间增量为我的信号“优化”滤波器。因此,对于上述信号,我可能希望在前 5 秒内过滤 580-620Hz 左右,然后在接下来的 5 秒内过滤 630-670Hz。本质上,我想最大限度地减少最终 IF 信号中的噪声。

根据您发送的示例代码,我编写了一个函数,该函数使用 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))

使用此过滤器与希尔伯特变换例程结合获得的 IF 数据与使用 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 与 scipy.signal 进行比较

我现在的问题是,我可以将 audiolazy Butterworth 例程与您在原始帖子中描述的随时间变化的属性结合起来吗?

4

3 回答 3

7

AudioLazy与时变过滤器一起工作

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 变换滤波器方程中直接应用时变系数。

编辑:有关 AudioLazy 组件的更多信息

以下示例是使用 IPython 完成的。

Resonator 是一个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)

其中Hz只是将单位从 Hz 更改为 rad/sample 的数字,如大多数 AudioLazy 组件中使用的那样。

让我们使用freq = pi/4and bw = pi/8(pi已经在audiolazy命名空间中):

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()Stream:

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 一维数组也是可迭代的:

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]

最后一个示例显示了时变行为。

编辑 3:分块重复序列行为

line 函数的行为类似于numpy.linspace,它获取范围“length”而不是“step”。

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]

并在第一个示例的行中使用它,因此freqandbw变为:

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)相反,这些应该做同样的事情。

编辑 5:时变巴特沃斯滤波器

我花了最后几个小时试图让巴特沃斯作为你的例子进行个性化,order = 2直接施加并给出半功率带宽(~3dB)。我已经做了四个示例,代码在这个 Gist中,并且我更新了 AudioLazy 以包含gauss_noise高斯分布的噪声流。请注意,gist 中的代码没有进行任何优化,仅在这种特殊情况下工作,并且由于“每个样本”系数查找行为,啁啾示例使其非常慢。瞬时频率可以从 rad/sample 中的 [过滤] 数据中获得:

diff(unwrap(phase(hilbert(filtered_data))))

在哪里diff = 1 - z ** -1或另一种方法在离散时间内找到导数,hilbert是函数来自scipy.signal它给我们分析信号(离散希尔伯特变换是其结果的虚部),另外两个是来自 AudioLazy 的辅助函数。

这就是当巴特沃斯在保持记忆的同时突然改变其系数而没有噪音时发生的情况:

variable_butterworth_abrupt_pure_sinusoid.png

在这种转变中,一种明显的振荡行为是显而易见的。您可以使用移动中值在较低频率侧“平滑”,同时保持突然过渡,但这不适用于较高频率。好吧,这就是我们对完美正弦曲线的期望,但是如果有噪声(很多噪声,高斯的标准偏差等于正弦波幅度),它变成:

variable_butterworth_abrupt_noisy.png

然后我试着用啁啾做同样的事情,正是这样:

variable_butterworth_pure_sinusoid.png

当在最高频率使用较低带宽进行滤波时,这显示出一种奇怪的行为。加上加性噪声:

variable_butterworth_noisy.png

要点中的代码也是AudioIO().play最后一个嘈杂的啁啾。

编辑 6:时变谐振器滤波器

在同一个 Gist 中添加了一个使用谐振器而不是 Butterworth 的示例。它们在纯 Python 中并且没有经过优化,但比butter在啁啾期间调用每个样本执行得更快,并且更容易实现,因为所有resonator策略都接受 Stream 实例作为有效输入。以下是级联两个谐振器(即二阶滤波器)的图:

reson_2_abrupt_pure_sinusoid.png reson_2_abrupt_noisy.png reson_2_pure_sinusoid.png reson_2_noisy.png

对于三个谐振器的级联(即三阶滤波器)也是如此:

reson_3_abrupt_pure_sinusoid.png reson_3_abrupt_noisy.png reson_3_pure_sinusoid.png reson_3_noisy.png

这些谐振器在中心频率处具有等于 1 (0 dB) 的增益,并且即使在根本没有任何滤波的情况下,转换中“突然纯正弦”图的振荡模式也会发生。

于 2013-08-28T02:05:41.590 回答
1

如果您使用切片提取信号的一部分,那么您实际上是在使用矩形窗口对数据进行窗口化,由于突然的不连续性,该矩形窗口在边缘“响起”。解决此问题的一种方法是使用响铃较少的窗口,例如汉宁窗口:

import numpy as np

signal = np.random.randn(222)
length = 50
window = np.hanning(length)
for i in range(0, len(signal)-length, length):
    do_filtering(signal[i:i+length] * window)

更多关于窗口:http ://en.m.wikipedia.org/wiki/Window_function

于 2013-08-08T15:44:29.303 回答
1

尝试使用您正在使用的每个滤波器对整个信号进行滤波,然后适当地合并滤波后的信号。粗略地说,在伪代码中:

# 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)。这是一个使用 numpy 的实现

基本上,您可以对信号进行 STFT,通过对时频阵列进行操作来过滤信号,然后对阵列进行逆 STFT 以获得过滤后的信号。

使用可逆小波变换可能会得到更好的结果。这是一篇付费论文,描述了如何使用小波变换完成几乎所有您尝试做的事情。

于 2013-08-08T14:16:16.920 回答