Os sinais do mundo real e Sistemas de Caso: Projeto filtro analógico com uma torção

Você é dado a tarefa de projetar uma entrada analógica (de tempo contínuo) filtrar para atender as especificações de resposta amplitude mostrados. Você também precisa encontrar a resposta ao degrau do filtro, determinar o valor da superação de pico e tempo onde ocorre a superação de pico.

O objectivo do desenho de filtro é para a magnitude da resposta de frequência em dB (20log₁₀|H(f) |) Para passar através da região não sombreada da figura como a frequência aumenta. Os requisitos de concepção reduzir a banda passante e stopband frequências críticas f_p e f_s Hz e os níveis de atenuação de banda passante e stopband UMA_p e UMA_s dB.

Além disso, a característica de resposta é para ser Butterworth, o que significa que a função de resposta e sistema de magnitude filtro de tomar esta forma:

Aqui, N é a ordem do filtro, f_c representa a banda de passagem de 3 dB frequência de corte do filtro, e os pólos, localizados em um semicírculo é a esquerda; meia s-plano, são dadas pela

Este problema exige que você trabalhe no domínio da frequência, o domínio do tempo, e talvez o s-domínio, dependendo da abordagem solução que você escolher.

A partir dos requisitos, a resposta de frequência do filtro tem ganho unitário (0 dB) na banda passante. A resposta de passo (a caracterização no domínio do tempo) do filtro Butterworth é conhecido por ultrapassar a unidade antes de finalmente se decidir a unidade como

Para projetar o filtro, você pode usar uma das duas abordagens:

1. Trabalhar uma solução à mão, usando a resposta de frequência magnitude Butterworth |H_BU(f) | e a função do sistema, H_BU(s).

2. Use os recursos de projeto do filtro do SciPy sinal pacote.

Encontrar a ordem do filtro e 3 dB frequência de corte

Siga estes passos para projetar o filtro usando Python e SciPy para fazer o número real de trituração:

1. Encontrar N e f_c para satisfazer os requisitos de resposta de magnitude.

   Use a função SciPy N, wc = signal.buttord (wp, ws, Ap, como, analógico = 1) e entre os requisitos de projeto do filtro, onde wp e ws são os de banda passante e stopband frequências críticas em rad / s e Ap e Como são os níveis de banda de passagem da faixa de rejeição e de atenuação (ambos os conjuntos de números vir a partir da figura precedente). A função retorna a ordem do filtro N e a frequência de corte banheiro em rad / s.

2. Sintetizar o filtro - encontrar o {b_k} E {uma_k} coeficientes da equação diferencial LCC que realiza o sistema desejado.

   Se encontrar elementos do circuito é o fim do jogo, você pode ir lá imediatamente, usando fórmulas de síntese circuito. Chamar a função SciPy b, a = signal.butter (N, wc, analógico = 1) com a ordem do filtro e a freqüência de corte, e ele retorna os coeficientes do filtro em matrizes b e uma.

3. Encontre a resposta ao degrau em forma matemática exata ou através de simulação.

Veja como utilizar as ferramentas de Python com as exigências do projeto dadas e, em seguida, verificar o trabalho traçando a resposta de frequência como uma sobreposição. Nota: Você pode fazer a mesma coisa em MATLAB com quase a mesma sintaxe.

Dentro [379]: N, wc = signal.buttord (2 * pi * * pi * 1e3,2 10e3,3.0,50, analógico = 1) ordem do filtro # find NIn [380]: N # filtro orderOut [380]: 3In [381]: Wc freq # corte em rad / sOut [381]: 9222.4701630595955In [382]: B, a = signal.butter (N, wc, analógico = 1) # get coeffs.In [383]: ataque[383]: Array ([7.84407571e + 11 + 0.j]) Em [384]: Real (a) Out [384]: Array ([1.00000000e + 00, 1.84449403e + 04,1.70107912e + 08,7.84407571e + 11])

Os resultados da Linha [379] dizer-lhe que a ordem do filtro necessário é N = 3 e a frequência de corte do filtro necessário é

Os conjuntos de coeficientes de filtro, também são incluídos nos resultados.

Use o real() função para exibir de forma segura a parte real da matriz de coeficientes uma porque você sabe que os coeficientes são reais. Como? Os pólos, raízes denominador H_BU(s), São reais ou ocorrem em pares conjugados complexos, garantindo que o polinômio denominador tem coeficientes reais quando multiplicado fora. As peças muito pequenas imaginários são devido a erros de precisão numérica.

Verificar a resposta final frequência projeto

Para verificar o design, usar a receita resposta de frequência.

Dentro [386]: F = logspace (2,5,500) axisIn # log frequência [387]: W, H = signal.freqs (b, a, 2 * pi * f) Em [388]: Semilogx (f, 20 * log10 (abs (H)), 'g')

A figura mostra o gráfico da resposta final design de magnitude, juntamente com os requisitos de design originais.

Encontrar a resposta ao degrau dos coeficientes do filtro

A abordagem mais elegante para encontrar a resposta ao degrau dos coeficientes do filtro é encontrar

o s-seção de domínio da figura diz-lhe como completar a expansão em frações parciais (PFE) numericamente. Você tem as matrizes de coeficiente para H(s), Então tudo que você precisa fazer é multiplicar o polinômio denominador por s. Você pode fazer isso manualmente ou você pode usar uma relação entre coeficientes polinomiais e convolução seqüência.

Quando você multiplicar dois polinômios, as matrizes de coeficientes para cada polinômio são convolved, como em convolução seqüência.

Aqui, você trabalha com o problema, usando signal.convolve para realizar a multiplicação polinomial no denominador. Para convencê-lo de que isso realmente funciona, considere multiplicação dos dois polinômios seguintes:

(x² + x + 1) (x + 1) = x³ + 2x² + 2x + 1

Se você convolve os conjuntos de coeficientes [1, 1, 1] e [1, 1] como arrays em Python, você recebe essa saída:

Dentro [418]: Signal.convolve ([1,1,1], [1,1]) Out [418]: Array ([1, 2, 2, 1])

Isto concorda com o cálculo manual. Para encontrar o PFE, conecte as matrizes de coeficientes b e convolve (a, [1,0]) para dentro R, P, K = resíduo (b, a). Os coeficientes de [1, 0] corresponde ao s-polinomial de domínio s + 0.

Dentro [420]: R, P, K = signal.residue (b, signal.convolve ([1,0], um)) em [421]: R # (resíduos) arranhar pequena errorsOut numérica [421]: Array ([1.0000e + 00 + 2.3343e-16j, # resíduo 0, parte imag 0-1.0000e + 00 + 1.0695e-15j, # resíduo 1, parte imag 01.08935e-15 -5.7735e-01J, # resíduo 2, parte real 01.6081e-15 + 5.7735e-01J]) # resíduo 3, parte real 0cm [422]: P # (pólos) Out [422]: Array ([0.0000 + 0.0000e + 00j, # pólo 0-9222,4702 -1.5454e-12j, # pólo 1, imag parte 0-4.611,2351 -7.9869e + 03J, # pólo 2-4611,2351 + 7.9869e + 03J]) # pólo 3In [423]: K # (de divisão longa) Out [423]: Array ([0. + 0.j]) # adequada racional, por isso não termos

Você tem quatro pólos: dois reais e um par conjugado complexo - um pouco de confusão para o trabalho através, mas é factível. Consulte a par transformar

para calcular a transformada inversa para todos os quatro termos.

Para os pólos conjugadas, os resíduos são também conjugados. Essa propriedade sempre se mantém.

Você pode escrever a transformada inversa dos termos conjugado pólo como senos e co-senos, usando a fórmula de Euler eo cancelamento de partes imaginárias na frente do cosseno e partes reais na frente da sine:

Juntando tudo, você começa y_passo(t) = você(t) - e^-9,222.47tvocê(t) -2 X 0,5774e^-4,611.74tsin (7,986.89t)você(t). Tendo esta forma é bom, mas você ainda precisa encontrar o máximo de função para t > 0 e o local máxima. Para fazer isso, representar graficamente a função e observar o máximo.

Uma abordagem mais directa é a utilização de simulação via signal.lsim e a receita de domínio de tempo. A entrada do sistema é um passo, de modo que a saída de simulação será a resposta de passo. A partir da resposta ao degrau simulado, você pode calcular a superação de pico numericamente e vê-lo em uma trama. O código de linha de comando IPython é

Dentro [425]: T = arange (0,0.002,1e-6) # passo menos de menor tempo de Constantin [426]: T, ys, x_state signal.lsim = ((b, a), os (len (t)), T) [428]: Plot (t * 1e3, ys)

Usando o tempo de matriz t ea matriz resposta ao degrau ys, você pode usar o max () e encontrar() funções para concluir a tarefa:

Dentro [436]: Max (real (ys)) # real para limpar num. errorsOut [436]: 1,0814651457627822 # pico superação is8.14% In [437]: Encontrar (reais (YS) == max (reais (ys))) Out [437]: Array ([534]) # find pico de ser no índice 534In [439]: T [534] * 1e3 # tempo no índice 534 em msOut [439]: 0,5339