使用 Python lfilter 过滤信号答案

【问题标题】：Filtering signal with Python lfilter使用 Python lfilter 过滤信号
【发布时间】：2014-03-19 17:39:15
【问题描述】：

我是 Python 新手，在过滤信号时完全卡住了。这是代码：

import numpy as np
import matplotlib.pyplot as plt
from scipy import signal

fs=105e6
fin=70.1e6

N=np.arange(0,21e3,1)

# Create a input sin signal of 70.1 MHz sampled at 105 MHz
x_in=np.sin(2*np.pi*(fin/fs)*N)

# Define the "b" and "a" polynomials to create a CIC filter (R=8,M=2,N=6)
b=np.zeros(97)
b[[0,16,32,48,64,80,96]]=[1,-6,15,-20,15,-6,1]
a=np.zeros(7)
a[[0,1,2,3,4,5,6]]=[1,-6,15,-20,15,-6,1]

w,h=signal.freqz(b,a)
plt.plot(w/max(w),20*np.log10(abs(h)/np.nanmax(h)))
plt.title('CIC Filter Response')

output_nco_cic=signal.lfilter(b,a,x_in)

plt.figure()        
plt.plot(x_in)
plt.title('Input Signal')
plt.figure()        
plt.plot(output_nco_cic)
plt.title('Filtered Signal')

还有情节：

如您所见，虽然滤波器传递函数是正确的，但输出却不正确。而且我不明白为什么我的代码不起作用。我在 Matlab 中编写了相同的代码，输出看起来还不错。

谢谢你的帮助！

【问题讨论】：

没人？，一条线索？
这甚至没有在我的机器上运行，我得到了除以零

标签： python scipy filtering signal-processing

【解决方案1】：

我不觉得这不适用于 Python 令人困惑，但我确实觉得它适用于 Matlab 令人困惑。

CIC 过滤器不适用于浮点数。

更新：

有趣的是，至少对于我拥有的 scipy 版本，lfilter 不适用于整数数组——我得到一个 NotImplemented 错误。这是一个 CIC 过滤器的 numpy 版本，它的速度大约是我机器上纯 Python 实现的两倍：

# Implements an in-memory CIC decimator using numpy.

from math import log
from numpy import int32, int64, array

def cic_decimator(source, decimation_factor=32, order=5, ibits=1, obits=16):

    # Calculate the total number of bits used internally, and the output
    # shift and mask required.
    numbits = order * int(round(log(decimation_factor) / log(2))) + ibits
    outshift = numbits - obits
    outmask  = (1 << obits) - 1

    # If we need more than 64 bits, we can't do it...
    assert numbits <= 64

    # Create a numpy array with the source
    result = array(source, int64 if numbits > 32 else int32)

    # Do the integration stages
    for i in range(order):
        result.cumsum(out=result)

    # Decimate
    result = array(result[decimation_factor - 1 : : decimation_factor])

    # Do the comb stages.  Iterate backwards through the array,
    # because we use each value before we replace it.
    for i in range(order):
        result[len(result) - 1 : 0 : -1] -= result[len(result) - 2 : : -1]

    # Normalize the output
    result >>= outshift
    result &= outmask
    return result

【讨论】：

【解决方案2】：

代码很好，lfilter 在它创建的 float64 数组上工作正常。但是分母多项式“a”的所有根都在 z = 1，这使得滤波器“条件稳定”。由于数值不准确，它最终会发散。并且 70.1 MHz 的输入信号在通带之外，所以它不会在输出中显示太多。如果您将输入更改为 0.701 MHz 左右，并将信号缩短为 1000 个样本而不是 21000 个样本，您会看到它按原样工作。尝试这些更改，您会看到之后会发生什么：鳍=70.1e4 N=np.arange(0,2000,1) （为了摆脱除以零的抱怨，请在 log10 中添加 1.0e-12）

要正确执行 CIC，您需要一个正确处理条件稳定极点的实现。

【讨论】：