Рассчитываем DSP фильтры при помощи SciPy

23 октября 2023

SciPy — это библиотека для языка Python, содержащая различные методы для инженерных и научных расчетов. При помощи SciPy можно считать интегралы, решать задачи оптимизации, обрабатывать изображения, и не только. Особый интерес представляют методы для расчета DSP-фильтров. Давайте разберемся, как ими пользоваться.

Нам понадобятся следующие зависимости:

pip install matplotlib numpy scipy

Чтобы не возникло конфликтов с другими проектами, не лишено смысла завести отдельную песочницу при помощи virtualenv. SciPy тесно связан с Matplotlib и NumPy. Эти библиотеки мы уже проходили. Не будем на них задерживаться.

Воспользуемся следующим шаблоном скрипта:

#!/usr/bin/env python3

# firdesign.py
# Aleksander Alekseev, 2023
# https://eax.me/scipy-dsp-filter-design/

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

sample_rate = 44100.0        # Sample rate, Hz
nyq_rate = sample_rate / 2.0 # Nyquist rate

# ========================================
# Create a low-pass FIR filter
numtaps = 61
cutoff_hz = sample_rate/4
taps = signal.firwin(numtaps, cutoff_hz/nyq_rate, window = 'hamming')
# ========================================

for t in taps:
    print("{},".format(t))

# Plot the filter coefficients

plt.figure(1)
plt.plot(taps, 'o-', linewidth=2)
plt.title('Impulse response')
plt.grid(True)

# Plot the magnitude response

w, h = signal.freqz(taps)

plt.figure(2)
plt.plot((w/np.pi)*nyq_rate, 20*np.log10(np.abs(h)+(1.0e-10)),
         linewidth=2)
plt.xlabel('Frequency (Hz)')
plt.ylabel('Gain (dB)')
plt.ylim(-95, 5)
plt.title('Magnitude Response')
plt.grid(True)

# Plot the phase response
# see https://dsp.stackexchange.com/a/50564/36360

plt.figure(3)
h = h*np.exp(1j*w*(numtaps//2))
plt.plot((w/np.pi)*nyq_rate, 180*np.angle(h)/np.pi, linewidth=2)
plt.xlabel('Frequency (Hz)')
plt.ylabel('Phase (deg)')
plt.title('Phase Response')
plt.grid(True)

plt.show()

Выделенный фрагмент будем заменять, чтобы получались различные фильтры.

В приведенном примере рассчитывается FIR фильтр, представляющий собой ФНЧ. Количество коэффициентов определено по знакомой нам эвристике:

numtaps = Atten/(22*(f_stop/f_sample - f_pass/f_sample))
numtaps = 60/(22*(11_000/44100-9_000/44100)) = 60.14

Вот АЧХ фильтра:

АЧХ FIR фильтра

… а вот импульсная характеристика, она же коэффициенты фильтра:

Импульсная характеристика FIR фильтра

Примерно половина коэффициентов равна нулю. Это потому что фильтр имеет частоту среза sample_rate/4. Еще можно заметить, что коэффициенты фильтра симметричны. Благодаря этому факту можно сократить количество умножений еще в два раза.

Альтернативный код:

# The desired width of the transition from pass to stop,
# relative to the Nyquist rate
width = 2000/nyq_rate

# The desired attenuation in the stop band, in dB
ripple_db = 60.0

# Compute the order and Kaiser parameter for the FIR filter
numtaps, beta = signal.kaiserord(ripple_db, width)

# The cutoff frequency of the filter.
cutoff_hz = 10000

# Use firwin with a Kaiser window to create a lowpass FIR filter
taps = signal.firwin(numtaps, cutoff_hz/nyq_rate,
                     window=('kaiser', beta))

Здесь мы получаем фильтр, имеющий 81 коэффициент. Его АЧХ лучше, однако вычислительная сложность существенно выше.

Еще при помощи метода firwin() можно рассчитать ФВЧ:

numtaps = 51
cutoff_hz = 10000
taps = signal.firwin(numtaps, cutoff_hz/nyq_rate,
                     window = 'hamming', pass_zero=False)

… полосовой фильтр:

numtaps = 51
taps = signal.firwin(numtaps, [5000/nyq_rate, 15000/nyq_rate],
                     window = 'hamming', pass_zero=False)

… и режекторный фильтр:

numtaps = 51
taps = signal.firwin(numtaps, [5000/nyq_rate, 15000/nyq_rate],
                     window = 'hamming')

Метод remez() также предназначен для расчета FIR фильтров. Но в отличие от firwin(), он может генерировать дифференциаторы, а также преобразователи Гильберта. И те и другие играют важную роль в SDR.

Преобразователь Гильберта — это фазовый фильтр, или all-pass filter, вносящий во входные сигналы фазовый сдвиг 90°. Рассчитаем его при помощи remez():

numtaps = 51
window = signal.windows.kaiser(numtaps, beta=8)
taps = signal.remez(numtaps, [0, nyq_rate], [1],
                    type='hilbert', fs=sample_rate)
taps = taps * window

АЧХ фильтра:

АЧХ преобразователя Гильберта

… ФЧХ фильтра:

ФЧХ преобразователя Гильберта

… и его коэффициенты:

Коэффициенты преобразователя Гильберта

Заметьте, что половина коэффициентов снова оказалась равна нулю. Это общее свойство преобразователей Гильберта с нечетным количеством коэффициентов. Тот факт, что фильтр подавляет постоянные сигналы, тоже не случаен. Никто не знает, что такое фазовый сдвиг постоянного сигнала.

Выход фильтра имеет временную задержку. Чтобы относительный фазовый сдвиг входа и выхода был 90°, вход нужно задерживать на numtaps//2 сэмплов. Сделать это можно, дописав numtaps//2 нулей в начале последовательности сэмплов. Это еще одна причина, почему в преобразовании Гильберта используют нечетное количество коэффициентов. Когда оно четное, задержка фильтра составляет нецелое число сэмплов.

Аргументы метода remez() задают АЧХ фильтра. При желании можно получить ФВЧ, ДПФ, и так далее с фазовым сдвигом 90°. В данном примере во всей полосе [0, nyq_rate] фильтр обладает единичным усилением [1].

Тот же фильтр может быть рассчитан и без remez():

numtaps = 51
window = signal.windows.kaiser(numtaps, beta=8)
taps = [(1/(np.pi*n))*(1 - np.cos(np.pi*n))
         for n in range(1, (numtaps//2)+1)]

# an equivalent formula
# taps = [2*np.sin(np.pi*n/2)**2/(np.pi*n)
#          for n in range(1, (numtaps//2)+1)]

# +90 deg phase shift
taps = np.append(np.append(np.flip(taps), [0]), np.negative(taps))

# -90 deg phase shift
# taps = np.append(np.negative(np.flip(taps)), np.append([0], taps))

taps = taps * window

Занимательный факт. Радиолюбители-конструкторы вместо фазовых сдвигов 0° и 90° часто используют пару фильтров со сдвигами +45° и −45°. Данный подход вычислительно неэффективен. Если вам все же нужен фазовый сдвиг на 45°, его можно получить, сложив сигнал с ним же, сдвинутым на 90°. По тому же принципу получаются любые другие фазовые сдвиги.

Для расчета IIR фильтров SciPy предлагает два метода. Первый — это iirfilter():

#!/usr/bin/env python3

# iirdesign.py
# Aleksander Alekseev, 2023
# https://eax.me/scipy-dsp-filter-design/

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

sample_rate = 44100.0
nyq_rate = sample_rate / 2.0

# ========================================
# Create a Chebyshev band-pass IIR filter
b, a = signal.iirfilter(11, [5000, 15000], fs = sample_rate, rp=0.5,
                        btype='band', ftype='cheby1')
# ========================================

print("=== a ({}) ===".format(a.size))
for t in a:
  print("{},".format(t))

print("=== b ({}) ===".format(b.size))
for t in b:
  print("{},".format(t))

# Plot the impulse response

idx = np.arange(0, 256)
x = 1.0*(idx == 0)
ir = signal.lfilter(b, a, x)

plt.figure(1)
plt.plot(ir, 'o-', linewidth=2)
plt.title('Impulse response')
plt.grid(True)

# Plot the frequency response

w, h = signal.freqz(b, a)

plt.figure(2)
plt.title('Magnitude Response')
plt.plot((w/np.pi)*nyq_rate, 20 * np.log10(abs(h)), linewidth=2)
plt.ylabel('Amplitude [dB]')
plt.xlabel('Frequency')
plt.grid(True)
plt.ylim(-95, 5)

# Plot zeros and poles

zeros, poles, _ = signal.tf2zpk(b, a)
for m in np.abs(poles):
  if m >= 1.0:
    print("IIR FILTER IS UNSTABLE!")
    break

circle = patches.Circle(
           (0,0), radius=1, fill=False,
           color='black', ls='solid', alpha=0.3)

fig, ax = plt.subplots()
plt.title('Poles and Zeroes')
ax.add_patch(circle)
plt.plot(zeros.real, zeros.imag, 'o', markersize=5,
         alpha=0.7, label = "Zeros")
plt.plot(poles.real, poles.imag, 'x', markersize=5,
         alpha=0.7, label = "Poles")
plt.legend(loc = 'upper right')

plt.show()

Здесь мы рассчитываем фильтр Чебышева с полосой пропускания 5 кГц .. 15 кГц и рябью в полосе пропускания 0.5 dB:

АЧХ IIR фильтра

Проверить стабильность IIR фильтра можно двумя способами. Во-первых, подать единичный импульс и посмотреть на выход фильтра. Скрипт использует для этого метод lfilter():

Импульсная характеристика IIR фильтра

Фильтр не уходит в генерацию, а значит, стабилен.

Во-вторых, можно посмотреть, куда попадают нули и полюса. Метод iirfilter() возвращает передаточную функцию фильтра в области комплексных частот в формате коэффициентов полинома. Метод tf2zpk() переводит ее из формата коэффициентов полинома в формат корней многочлена, они же нули и полюса (исходный текст содержал ошибку, спасибо Kirill Kotyagin за уточнение):

Полюса и нули IIR фильтра

Чтобы фильтр был стабильным, полюса должны находиться внутри единичного круга на комплексной плоскости. Если полюс попал прямо на круг или вышел за его пределы, значит фильтр нестабилен. Нули могут находиться где угодно. Их рисуют просто так, для красоты. Видим, что фильтр стабилен. В качестве эксперимента допишите в скрипт a[0] *= 1.003 и сравните результат.

Второй метод для расчета IIR фильтров — это iirdesign():

# ...
b, a = signal.iirdesign(wp = 9000, ws = 11000, fs = sample_rate,
                        gpass = 1, gstop = 60)
# ...

Здесь рассчитывается ФНЧ с полосой пропускания 9 кГц, подавлением 60+ dB выше 11 кГц и потерями в полосе пропускания не более 1 dB. Данный метод удобен тем, что сам определяет требуемое количество коэффициентов. Их здесь получилось 16 штук. FIR фильтру для достижения схожей АЧХ требуется порядка 60 .. 80 коэффициентов. Данный IIR фильтр стабилен, как и предыдущий.

Таким образом, мы установили, что SciPy действительно имеет все необходимое для расчета DSP фильтров.

Метки: , , , .


Вы можете прислать свой комментарий мне на почту, или воспользоваться комментариями в Telegram-группе.

Копирование материалов данного сайта не возбраняется при условии
указания ссылки на первоисточник. © 2009–2024 Записки программиста