iPython notebook как альтернатива Matlab

[БендеровецЪ]

А вот пример как запилить коректирующий IIR по типу "линквица".

Для начала нада запилить фунцию конвертации аналогова фильтра в IIR. Там есть свои приколы, и даже несколько вариантов имплементации (импулс инвариант, билинейный).
Подробнее обсуждается сдесь: https://dsp.stackexchange.com/questions/...pole/19211
Примеры: http://www.musicdsp.org/files/Audio-EQ-Cookbook.txt

Ну и тут как всегда удобно сделать свою библиотечку (во вложении). Принимает три параметра: вектор b вектор a и 1/fs синтезируемого iir

Использование выглядит приблизительно вот так:

Код:

from analog2iir import analog_iir
b_z_cor, a_z_cor = analog_iir(a_s, a_s_new, 1./fs)

(для моего примера fs = 96000.)

Собственно предидущая строка уже содержит то как делается коректирующий фильтр для системы (b_s, a_s) для получения желаемого отклика аналогичному для системы с (b_s_new, a_s_new) прт одинаковых a.

Собсна если кто помнит как делается "коректор Линквица". В числитель ставится знаменатель имеющейся системы, а в знаменатель - знаменатель целевой. Что я и сделал сразу вписав "a_s" и "a_s_new" в соотв места.

Для проверки хорошо бы исследовать АЧХ и ФЧХ исходной системы, целефой системы и корректирующего фильтра. Удобно использоваю уже имеюшиеся функции freqs и freqz

Код:

w_s, h_s = sp.signal.freqs(b_s, a_s, worN=np.logspace(1, 5, 1000)*2.*np.pi)
w_s_new, h_s_new = sp.signal.freqs(b_s_new, a_s_new, worN=np.logspace(1, 5, 1000)*2.*np.pi)
w_z, h_z = signal.freqz(b_z_cor, a_z_cor, worN=np.logspace(1, 5, 1000)*2.*np.pi/fs)

строим все это безобразие на графике:

Код:

plt.figure()
plt.plot(w_s/(2*np.pi), np.absolute(h_s))
plt.plot(fs*w_z/(2*np.pi), np.absolute(h_z))
plt.plot(w_s/(2*np.pi), np.absolute(h_s_new))
plt.xscale('log')
plt.yscale('log')
plt.xlim(10, 1e5)
plt.ylim(1e-5, 1)
plt.grid(which = 'both')
plt.show()

plt.figure()
plt.plot(w_s/(2*np.pi), np.angle(h_s))
plt.plot(w_s/(2*np.pi), np.angle(h_s_new))
plt.plot(fs*w_z/(2*np.pi), np.angle(h_z))
plt.xscale('log')
plt.xlim(10, 1e5)
plt.grid(which = 'both')
plt.show()

и получаем



Для примера использовал следующие вводные:

Код:

fc = 500. #Hz
fc_new = 100. #Hz
w_s = 2.*np.pi*fc
w_s_new = 2.*np.pi*fc_new
q_s = 0.7
q_s_new = 0.7

b_s = np.array((np.square(1./w_s) , 0, 0))
a_s = np.array((np.square(1./w_s) , 1./(w_s*q_s) , 1))

b_s_new = np.array((np.square(1./w_s_new) , 0, 0))
a_s_new = np.array((np.square(1./w_s_new) , 1./(w_s_new*q_s_new) , 1))

При желании можно просимулировать сисемы во временной области (все эти ваши импульсные отклики или ступеньки) ипользуя готовые функции lsim и dlsim.

NB ну если где-то с чем-то промахнулся то просьба не пинать, это все-таки не моя специальность :)