python - pass - scipy help

Filtrando la señal con Python lfilter

Question

Soy nuevo con Python y estoy completamente atascado al filtrar una señal. Este es el código:

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

Y las tramas:

Como puede ver, aunque la función de transferencia de filtro es correcta, la salida no lo es. Y no puedo ver por qué mi código no funciona. He codificado lo mismo en Matlab y la salida se ve bien.

Thaks por la ayuda!

Answer 1

No me parece confuso que esto no funcione con Python, pero me resulta confuso que funcionó con Matlab.

Los filtros CIC no funcionan con números de coma flotante.

ACTUALIZAR:

Curiosamente, al menos con la versión de scipy que tengo, lfilter no funciona con matrices de enteros, obtengo un error NotImplemented. Aquí hay una versión numpy de un filtro CIC que es aproximadamente el doble de rápido que una implementación pura de Python en mi máquina:

# 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

Answer 2

El código está bien, y lfilter funciona bien en las matrices float64 que crea. Pero el denominador polinomial "a" tiene todas sus raíces en z = 1, lo que hace que el filtro sea "condicionalmente estable". Debido a la inexactitud numérica, eventualmente divergerá. Y la señal de entrada a 70.1 MHz está muy por fuera de la banda de paso, por lo que no aparece mucho en la salida. Si cambia la entrada a 0.701 MHz o menos, y acorta la señal a 1000 muestras en lugar de 21000, verá que funciona tal como está. Pruebe estos cambios y verá lo que sucede después: fin = 70.1e4 N = np.arange (0,2000,1) (y para deshacerse de la división por cero, agregue 1.0e-12 dentro del log10)

Para hacer un CIC correcto, necesita una implementación que maneje correctamente los polos condicionalmente estables.