- Educación
- Ciencia
- Ingeniería
- Caso de Señales y Sistemas del Mundo Real: Diseño de filtro analógico con un giro
Libro Relacionado
TEMARIO
Señales y sistemas para Explicado
Por Mark Wickert
Se le asigna la tarea de diseñar un filtro analógico (de tiempo continuo) para cumplir con las especificaciones de respuesta de amplitud que se muestran. También es necesario encontrar la respuesta del paso del filtro, determinar el valor de la sobreoscilación del pico y el momento en que se produce la sobreoscilación del pico.
El objetivo del diseño del filtro es que la magnitud de la respuesta de frecuencia en dB (20log10|H(f)|) pase a través de la región no sombreada de la figura a medida que aumenta la frecuencia. Los requisitos de diseño reducen a la banda de paso y a la banda de parada las frecuencias críticas fp y fs Hz y los niveles de atenuación de la banda de paso y de parada Ap y As dB.
Además, la característica de respuesta debe ser Butterworth, lo que significa que la respuesta de magnitud de filtro y la función del sistema toman esta forma:
Aquí, N es el orden del filtro, fc es la frecuencia de corte de la banda de paso de 3 dB del filtro, y los polos, situados en un semicírculo es la mitad izquierda del plano s, están dados por
Este problema requiere que trabaje en el dominio de frecuencia, el dominio de tiempo y quizás el s-dominio, dependiendo del enfoque de solución que elija.
A partir de los requisitos, la respuesta de frecuencia del filtro tiene ganancia unitaria (0 dB) en la banda de paso. La respuesta escalonada (una caracterización de tiempo-dominio) del filtro Butterworth es conocida por sobrepasar la unidad antes de que finalmente se asiente a la unidad como
Para diseñar el filtro, puede utilizar uno de dos enfoques:
- Trabaje una solución a mano, utilizando la respuesta de frecuencia de magnitud de Butterworth |HBU(f)| y la función de sistema, HBU(s).
- Utilice las capacidades de diseño de filtros del paquete de señales de SciPy.
Búsqueda del orden del filtro y de la frecuencia de corte de 3 dB
Siga estos pasos para diseñar el filtro usando Python y SciPy para hacer el cálculo numérico real:
- Use la función SciPy N,wc=signal.buttord(wp,ws,Ap,As,analog=1) e introduzca los requisitos de diseño del filtro, donde wp y ws son las frecuencias críticas de paso y banda de parada en rad/s y Ap y As son los niveles de atenuación de paso y banda de parada (ambos conjuntos de números provienen de la figura anterior). La función devuelve el orden del filtro N y la frecuencia de corte wc en rad/s.
- Sintetizar el filtro – encontrar los coeficientes {bk} y {ak} de la ecuación diferencial de LCC que realiza el sistema deseado Si encontrar elementos de circuito es el juego final, puede ir allí inmediatamente, utilizando fórmulas de síntesis de circuito. Llamar a la función SciPy b,a=signal.butter(N,wc,analog=1) con el orden del filtro y la frecuencia de corte, y devuelve los coeficientes de filtro en los arrays b y a.
- Encuentre la respuesta al paso en forma matemática exacta o mediante simulación.
He aquí cómo utilizar las herramientas Python con los requisitos de diseño dados y luego comprobar el trabajo trazando la respuesta de frecuencia como una superposición. Nota: Puede hacer lo mismo en MATLAB con casi la misma sintaxis.
En[379]: N,wc =signal.buttord(2*pi*1e3,2*pi*10e3,3.0, 50,analog=1) # buscar filtro orden NIn[380]: N # filter orderOut[380]: 3In[381]: wc # frecuencia de corte en rad/sOut[381]: 9222.470163059595955In[382]: b,a = signal.butter(N,wc,analog=1) # get coeffs.En[383]: bOut[383]: array ([7.84407571e+11+0.j])En[384]: real(a)Out[384]: array([1.00000000e+00, 1.84449403e+04, 1.70107912e+08,7.84407571e+11])
Los resultados de la Línea[379] le indican que el orden de filtrado requerido es N = 3 y la frecuencia de corte de filtro requerida es
Los conjuntos de coeficientes de filtro también se incluyen en los resultados.
Utilice la función real() para mostrar con seguridad la parte real de la matriz de coeficientes a porque sabe que los coeficientes son reales. ¿Cómo? Los polos, raíces del denominador de las UHB, son reales u ocurren en pares conjugados complejos, asegurando que el polinomio del denominador tenga coeficientes reales cuando se multiplica. Las partes imaginarias muy pequeñas se deben a errores de precisión numérica.
Comprobación de la respuesta en frecuencia final de diseño
Para comprobar el diseño, utilice la receta de respuesta de frecuencia.
En[386]: f = logspace(2,5,500) # log frequency axisIn[387]: w,H = signal.freqs(b,a,2*pi*f)En[388]: semilogx(f,20*log10(abs(H)),'g')
La figura muestra el diagrama de la respuesta de la magnitud de diseño final junto con los requisitos de diseño originales.
Encontrar la respuesta de paso a partir de los coeficientes de filtro
El enfoque más elegante para encontrar la respuesta de paso a partir de los coeficientes de filtro es encontrar
La sección s-domain de la figura le indica cómo completar la expansión de la fracción parcial (PFE) numéricamente. Tienes las matrices de coeficientes para H(s), así que todo lo que necesitas hacer es multiplicar el polinomio denominador por s. Puedes hacer esto a mano o puedes usar una relación entre los coeficientes polinomiales y la convolución de la secuencia.
Cuando se multiplican dos polinomios, los arrays de coeficiente para cada polinomio son convolucionados, como en la convolución secuencial.
Aquí se trabaja el problema, usando signal.convolve para realizar la multiplicación polinómica en el denominador. Para convencerte de que esto realmente funciona, considera la multiplicación de los dos polinomios siguientes:
Si se convuelven los coeficientes[1, 1, 1, 1] y[1, 1] como matrices en Python, se obtiene esta salida:
En[418]: señal.convolve([1,1,1],[1,1])Out[418]: array([1, 2, 2, 1])
Esto concuerda con el cálculo de la mano. Para encontrar la ZFP, enchufe los coeficientes de las matrices b y convolve(a,[1,0]) en R,p,K = residuo(b,a). Los coeficientes[1, 0] corresponden al polinomio s-dominio s + 0.
En[420]: R,p,K = residuo de señal (b,convolución de señal([1,0],a))En[421]: R #(residues) scratch tiny numerical errorsOut[421]:array([ 1.0000e+00 +2.3343e-16j, # residue 0, imag part 0 -1.0000e+00 +1.0695e-15j, # residue 1, imag part 0 1.08935e-15 -5.7735e-01j, # residue 2, real part 0 1.6081e-15 +5.7735e-01j]) # residuo 3, parte real 0In[422]: p #(polos)Out[422]:array([ 0.0000 +0.0000e+00j, # polo 0 -9222.4702 -1.5454e-12j, # polo 1, parte imag 0 -4611.2351 -7.9869e+03j, # polo 2 -4611.2351 +7.9869e+03j])# polo 3In[423]: K #(de la división larga)Out[423]: array([ 0.+0.j]) # racional adecuado, así que sin términos
Tienes cuatro polos: dos pares de conjugados reales y uno complejo, un poco desordenado, pero es factible. Consulte el par de transformación
para calcular la transformación inversa de los cuatro términos.
Para los polos conjugados, los residuos son también conjugados. Esta propiedad siempre se mantiene.
Se puede escribir la transformación inversa de los términos del polo conjugado como senos y cosenos, usando la fórmula de Euler y la cancelación de las partes imaginarias frente al coseno y las partes reales frente al seno:
Poniendo todo junto, usted obtiene ystep(t) = u(t) – e-9,222.47tu(t)-2 x 0.5774e-4,611.74tsin(7,986.89t)u(t). Tener este formulario es bueno, pero aún necesita encontrar el máximo de la función para t > 0 y la ubicación máxima. Para ello, dibuje la función y observe el máximo.
Un enfoque más directo es utilizar la simulación a través de signal.lsim y la receta del dominio de tiempo. La entrada del sistema es un paso, por lo que la salida de la simulación será la respuesta al paso. A partir de la respuesta de paso simulada, se puede calcular numéricamente el rebasamiento del pico y verlo en un gráfico. El código de la línea de comandos de IPython es
En[425]: t = arange(0,0.002,1e-6) # step less than smallest time constantIn[426]: t,ys,x_state = signal.lsim((b,a),ones(len(t)),t)En[428]: plot(t*1e3,ys)
Usando la matriz de tiempo t y la matriz de respuesta de pasos ys, puede usar las funciones max() y find() para completar la tarea:
En[436]: max(real(ys)) # real to clear num. errorsOut[436]: 1.0814651457627822 # sobregiro de pico es8.14%In[437]: find(real(ys)== max(real(ys)))Out[437]: array([534]) # find peak to be at index 534In[439]: t[534]*1e3 # time at index 534 in msOut[439]: 0.5339