2014-02-19 33 views
5

我是Python新手,過濾信號時完全卡住了。這是代碼:使用Python lfilter過濾信號

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

而且情節:

Input, filter and output signals

正如你可以看到,雖然濾波器的傳遞函數是正確的,輸出不是。而我不明白爲什麼我的代碼無法正常工作。我在Matlab中編碼相同,輸出看起來不錯。

Thaks的幫助!

+0

沒有人?,一個線索? – Alex

+0

這甚至沒有在我的機器上運行,我得到了除零 – Nate

回答

6

我不覺得困惑,這不適用於Python,但我確實發現它與Matlab一起工作令人困惑。

CIC過濾器不適用於浮點數。

UPDATE:

有趣的是,至少有SciPy的版本我都有,lfilter不整數數組工作 - 我得到一個NotImplemented錯誤。這裏是一個CIC濾波器也差不多快兩倍,我的機器上的純Python實現的numpy的版本:

# 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 
1

的代碼是好的,lfilter會在它創建float64數組罰款。但分母多項式「a」的所有根都在z = 1,這使得濾波器「條件穩定」。由於數字不準確,它最終會發生分歧。而70.1 MHz的輸入信號遠離通帶,因此在輸出中不會顯示太多。如果將輸入更改爲0.701 MHz或其附近,並將信號縮短爲1000個採樣而不是21000,則會看到它的工作原理。嘗試這些更改,你會看到之後會發生什麼: fin = 70.1e4 N = np.arange(0,2000,1) (並且爲了擺脫零投訴的分歧,添加1.0e-12內部log10)

要做CIC的權利,你需要一個實現,正確地處理條件穩定的極點。