Установка Wn для аналоговой модели Бесселя в scipy.signal

Question

Я пытался реализовать аналоговый фильтр Бесселя с частотой среза 2 кГц, используя scipy.signal, и меня смущает то, какое значение Wn нужно установить, поскольку в документации для состояний Wn (для аналоговых фильтров) должно быть установлено значение * Частота 1021 * (приблизительно 12000 рад / с). Но если я добавлю это к моим 1-секундным фиктивным данным с полусекундным импульсом, дискретизированным с частотой 500 000 Гц, я получу строку из 0 с и нанов. Чего мне не хватает?

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

def make_signal(pulse_length, rate = 500000):
    new_x = np.zeros(rate)
    end_signal = 250000+pulse_length
    new_x[250000:end_signal] = 1
    data = new_x
    print (np.shape(data))
    # pad on both sides
    data=np.concatenate((np.zeros(rate),data,np.zeros(rate)))
    return data 

def conv_time(t):
    pulse_length = t * 500000
    pulse_length = int(pulse_length)
    return pulse_length

def make_data(ti): #give time in seconds
    pulse_length=conv_time(ti)
    print (pulse_length)
    data = make_signal(pulse_length)
    return data
time_scale = np.linspace(0,1,500000)
data = make_data(0.5)    


[b,a] = scipy.signal.bessel(4, 12566.37, btype='low', analog=True, output='ba', norm='phase', fs=None)

output_signal = scipy.signal.filtfilt(b, a, data)
plt.plot(data[600000:800000])

plt.plot(output_signal[600000:800000])

При построении ответа с использованием freqs это не кажется мне плохим; где я ошибаюсь?

Answer 1

Вы передаете аналоговый фильтр функции scipy.signal.filtfilt, которая ожидает цифровой (т. Е. Дискретное время) фильтр. Если вы собираетесь использовать filtfilt или lfilter, фильтр должен быть цифровым.

Для работы с системами с непрерывным временем взгляните на функции

• scipy.signal.impulse (scipy.signal.impulse2)
• scipy.signal.step (scipy.signal.step2)
• scipy.signal.lsim (scipy.signal.lsim2)

(версии 2 решают ту же математическую задачу, что и версия без 2, но используют другой метод. В большинстве случаев версия без 2 подходит и на намного быстрее, чем 2 версия.)

Другие связанные функции и классы перечислены в разделе Линейные системы с непрерывным временем документации SciPy.

Например, вот сценарий, который отображает импульсные и пошаговые ответы вашего фильтра Бесселя:

import numpy as np
from scipy.signal import bessel, step, impulse
import matplotlib.pyplot as plt


order = 4
Wn = 2*np.pi * 2000
b, a = bessel(order, Wn, btype='low', analog=True, output='ba', norm='phase')

# Note: the upper limit for t was chosen after some experimentation.
# If you don't give a T argument to impulse or step, it will choose a
# a "pretty good" time span.
t = np.linspace(0, 0.00125, 2500, endpoint=False)
timp, yimp = impulse((b, a), T=t)
tstep, ystep = step((b, a), T=t)


plt.subplot(2, 1, 1)
plt.plot(timp, yimp, label='impulse response')
plt.legend(loc='upper right', framealpha=1, shadow=True)
plt.grid(alpha=0.25)
plt.title('Impulse and step response of the Bessel filter')

plt.subplot(2, 1, 2)
plt.plot(tstep, ystep, label='step response')
plt.legend(loc='lower right', framealpha=1, shadow=True)
plt.grid(alpha=0.25)
plt.xlabel('t')
plt.show()

Сценарий генерирует этот график: