Ответ IIR в Python

У меня есть новый вопрос, напрямую связанный с этой записью. – Встроенный в Python. У меня есть полосовой фильтр второго порядка IIR с заданными характеристиками. [Следующий код преднамеренно идиоматичен]:

fs = 40e6 # 40 MHz f sample frequency fc = 1e6/fs # 1 MHz center cutoff BW = 20e3/fs # 20 kHz bandwidth fl = (fc - BW/2)/fs # 0.99 MHz f lower cutoff fh = (fc + BW/2)/fs # 1.01 MHz f higher cutoff 

Который дает коэффициенты:

 R = 1 - (3*BW) K = (1 - 2*R*np.cos(2*np.pi*fc) + (R*R)) / (2 - 2*np.cos(2*np.pi*fc)) a0 = 1 - K # a0 = 0.00140 a1 = 2*(KR)*np.cos(2*np.pi*fc) # a1 = 0.00018 a2 = (R*R) - K # a2 = -0.00158 b1 = 2*R*np.cos(2*np.pi*fc) # b1 = 1.97241 b2 = -(R*R) # b2 = -0.99700 

Как было предложено ukrutt в предыдущем посте, я использовал scipy.signal.freqz, но, к сожалению, не получил ответа, который я искал, – который сказал, что я считаю, что фильтр работает по назначению (код ниже). Вот результат freqz:

введите описание изображения здесь

Мой вопрос: как я могу сгенерировать график больше, чем предполагаемый ответ?

Код:

 a = [0.0014086232031758072, 0.00018050359364826498, -0.001589126796824103] b = [1.9724136161684902, -0.9970022500000001] w,h = signal.freqz(a, b) h_dB = 20 * np.log10(np.abs(h)) plt.plot(w/np.max(w),h_dB) plt.grid() 

    Вы не увидите ничего симпатичного с линейной шкалой x. Я не знаю numpy, но я знаком с Matlab, и есть некоторые функции, которые можно делать в журналах. Попробуйте использовать масштаб x-log с:

     import matplotlib.pyplot as pyplot fig = pyplot.figure() ax = fig.add_subplot(2,1,1) line, = ax.plot(w/np.max(w), h_dB, color='blue', lw=2) ax.set_xscale('log') show() 

    Я не тестировал его, я не установил python 🙁

    Редактировать:

    Я попытался смоделировать фильтр butterworth в matlab для порядка 4 фильтра IIR и одного порядка 20 фильтра IIR.

     %!/usr/local/bin/matlab %% Inputs fs = 40e6; fc = 1e6; BW = 20e3; fl = (fc - BW/2); fh = (fc + BW/2); %% Build bandpass filter IIR Butterworth order 4 N = 4; % Filter Order h = fdesign.bandpass('N,F3dB1,F3dB2', N, fl, fh, fs); Hd1 = design(h, 'butter'); %% Build bandpass filter IIR Butterworth order 50 N = 20; % Filter Order h = fdesign.bandpass('N,F3dB1,F3dB2', N, fl, fh, fs); Hd2 = design(h, 'butter'); %% Compare fvtool(Hd1,Hd2); 

    Два фильтра

    Небольшой зум

    И здесь коэффициенты A и B для первого фильтра:

     FilterStructure: 'Direct-Form II Transposed' A: [2.46193004641106e-06 0 -4.92386009282212e-06 0 2.46193004641106e-06] B: [1 -3.94637005453608 5.88902106889851 -3.93761314372475 0.995566972065978] 

    Если я получу какое-то время, я попытаюсь сделать то же самое с numpy!

    Я не думаю, что проблема заключается в том, как вы анализируете ответ – это в вашем выборе фильтра. Вы пытаетесь создать очень узкий ответ фильтра, используя только фильтр низкого порядка IIR. Я думаю, вам нужен либо фильтр более высокого порядка, либо ослаблять ваши ограничения.

    Например, следующее использует фильтр butterworth, реализованный как IIR, который дает ответ, который имеет форму, более похожую на то, что вы ищете. Очевидно, что для получения ожидаемой характеристики фильтра потребуется больше работы.

      b, a = signal.butter(4, [1.0/4-1.0/2e2,1.0/4+1.0/2e2], 'bandpass', analog=False) w, h = signal.freqs(b, a) import matplotlib.pyplot as plt fig = plt.figure() plt.title('Digital filter frequency response') ax1 = fig.add_subplot(111) plt.semilogy(w, np.abs(h), 'b') plt.ylabel('Amplitude (dB)', color='b') plt.xlabel('Frequency (rad/sample)') ax2 = ax1.twinx() angles = np.unwrap(np.angle(h)) plt.plot(w, angles, 'g') plt.ylabel('Angle (radians)', color='g') plt.grid() plt.axis('tight') plt.show() 

    который дает:

    Участок реакции фильтра

    Проблема в том, что signal.freqz возвращает точки на полукруге … поэтому вы не можете расширять до большего диапазона x , если только вы не делаете это через signal.freqz . Я попытался немного подтолкнуть его, и я увидел, что вы можете использовать pass whole=True для signal.freqz , и вы получите то, что у вас есть выше, но зеркально отразились на отрицательном x . Вот и все. Тем не менее, есть еще один аргумент ключевого слова, который позволяет вам передать массив x точек, которые вы хотите signal.freqz на … поэтому я попытался использовать np.arange(-5., 5., 0.1) … и он не выглядел на всех, как сюжет, который вы ожидаете справа, это выглядело как куча отражений исходного сюжета. Так что это заставляет меня думать … Может быть, у сюжета, который у вас есть справа, а у левой – разные топоры? В частности, это одна угловая частота, а другая простая старая частота?

    При дальнейшем signal.freqz возвращает w,h , где w – нормализованная угловая частота в радианах / образце. Таким образом, вам не нужно выполнять нормализацию по np.max(w) в вашем коде, чтобы сделать сюжет. Тем не менее это все еще не решает проблему. Ваш график справа находится в единицах fc , а fc – в МГц (например, 1 / образец).

    Поэтому, чтобы сделать график слева согласным справа, я предполагаю, что это означает, что вам нужно «ненормализовать» вашу ось x , тогда вам нужно будет преобразовать данные из единиц угловой частоты, используемой в МГц.

    Или, может быть, более вероятно, используйте другую функцию, чем signal.freqz .