Давно ничего не писал и, может быть, так бы продолжалось и дальше, если бы не комментарий читателя с ником «Pashka41″ к статье под названием «Цифровой БИХ фильтр (IIR filter). Пример программы.». Статья, в порыве страсти, была написана мною довольно таки давно и содержит реализацию БИХ фильтров разнобразных видов. Комментарий пользователя содержал вопрос относительно типа фильтра использованного в реализации:

… Подскажите, какой тип фильтра вы выбрали? Чебышев, Баттерворт..? Какого порядка?

Прежде чем ответить на комментарий, я поставил себя на место «Pashka41″. Дело в том, что формулы для реализации были взяты мною из просторов «великого и могучего» и вспомнить откуда именно, удалось не сразу. Поэтому, закатав рукава, я начал смотреть мою С++ реализацию, делать обратное Z преобразование, обратное Билинейное преобразование (зачем-то) и тут здравый смысл (некоторые называют это ленью) вернулся ко мне.

Нет, вы не подумайте. Если бы я был жителем необитаемого острова и мне позарец нужно было бы найти ответ на данный вопрос, то я бы его искал дальше и не сдавался. К счастью, я всего-лишь обитатель космического корабля под названием «Земля». Поэтому, пошёл иным путём.

Я открыл поисковик и нашёл материал, который я использовал. К счастью, этот материал содержит всю необходимую мне (и читателю) информацию. Использованный мною в статье фильтр называют «Биквадратным» фильтром (eng. biquadratic filter) второго порядка.

Мне стало ужасно неудобно (отчасти даже стыдно) за такое халатное отношение к использованному материалу и я решил хоть как-то реабилитироваться. Поскольку в приводимой статье присутствуют передаточные функции для каждого из фильтров, я решил исследовать их АЧХ и ФЧХ, а результатами поделиться с вами. В конце концов набросал простенький скрипт на Python

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
#!/usr/bin/env python3
# -*- coding: utf-8 -*- 
""" 
Created on Wed May 16 21:52:15 2018 
@author:  Viktor Signaievskyi 
@blog:    www.HowToDoIt.com.ua 
@e-mail:  piligrim2007@meta.ua 
@details: Program which draws frequency responses for all kind of filters taken from
          http://www.musicdsp.org/files/Audio-EQ-Cookbook.txt
Copyrighted by ©Viktor Signaievskyi 
"""
 
import numpy as np 
import matplotlib.pyplot as plot 
from enum import Enum
 
omega = np.arange(0, np.pi, 0.01);
s = 1j*omega;
 
class FILTER_TYPE(Enum):
    LPF = 1
    HPF = 2
    BPF1 = 3
    BPF2 = 4
    NOTCH = 5
    APF = 6
    PEAKING_EQ = 7
    LOW_SHELF = 8
    HIGH_SHELF = 9
 
def get_filter_transfer_function(type):
    K = 0
    A = 0.1;
    Q = 1.0/np.sqrt(2.0);
    
    if type == FILTER_TYPE.LPF:
        K = 1.0 / (s*s + s/Q + 1.0);
    
    elif type == FILTER_TYPE.HPF:
        K = (s*s) / (s*s + s/Q + 1.0);
    
    elif type == FILTER_TYPE.BPF1:
        K = (s) / (s*s + s/Q + 1.0);
    
    elif type == FILTER_TYPE.BPF2:
        K = (s/Q) / (s*s + s/Q + 1.0);
    
    elif type == FILTER_TYPE.NOTCH:
        K = (s*s + 1.0) / (s*s + s/Q + 1.0);
    
    elif type == FILTER_TYPE.APF:
        K = (s*s - s/Q + 1.0) / (s*s + s/Q + 1.0)
    
    elif type == FILTER_TYPE.PEAKING_EQ:
        
        K = (s*s + s*(A/Q) + 1) / (s*s + s/(A*Q) + 1.0)
    
    elif type == FILTER_TYPE.LOW_SHELF:
        K = A * (s*s + (np.sqrt(A)/Q)*s + A)/(A*s*s + (np.sqrt(A)/Q)*s + 1.0);
    
    elif type == FILTER_TYPE.HIGH_SHELF:
        K = A * (A*s*s + (np.sqrt(A)/Q)*s + 1.0)/(s*s + (np.sqrt(A)/Q)*s + A);
        
    return K;
 
def draw_frequency_response(type):
 
    K = get_filter_transfer_function(type);
    
    if type == FILTER_TYPE.APF:
        axes = plot.gca()
        axes.set_ylim([-1, 1])
 
    plot.plot(omega, 20*np.log10(np.abs(K)[:]), "b-");
    plot.xscale('log')
    plot.xlabel("Frequency, Hz")
    plot.ylabel("|K(jw)|, dB")
    plot.grid()
 
    return
 
def draw_phase_response(type):
  
    K = get_filter_transfer_function(type);
    
    plot.plot(omega, np.arctan(np.imag(K)[:]/np.real(K)[:]), "r-");
    plot.xscale('log')
    plot.xlabel("Frequency, Hz")
    plot.ylabel("Phase, rad")
    plot.grid()
    
    return
 
# drawing frequency and phase responses for all filter types
#for f in [FILTER_TYPE.APF]:
for f in FILTER_TYPE:
    filter_name = str(f).replace('FILTER_TYPE.', '');
    
    draw_frequency_response(f);
    plot.title(filter_name + ' (Frequency response)');
    plot.draw();
    plot.savefig(filter_name + '_afr', bbox_inches='tight', dpi=100, transparent=True);
    plot.show();
    
    draw_phase_response(f);
    plot.title(filter_name + ' (Phase response)');
    plot.draw();
    plot.savefig(filter_name + '_pfr', bbox_inches='tight', dpi=100, transparent=True);
    plot.show();

и получил следующие картинки:

LPF_afrLPF_pfr

HPF_afrHPF_pfr

BPF1_afr BPF1_pfr

 BPF2_afr BPF2_pfr

NOTCH_afr NOTCH_pfr

APF_afr APF_pfr

PEAKING_EQ_afr PEAKING_EQ_pfr

LOW_SHELF_afr LOW_SHELF_pfr

HIGH_SHELF_afr HIGH_SHELF_pfr

Теперь фильтрами можна пользоваться без боязни сделать что-то непоправимое :-) .

Ну, а если серьезно и без шуток, в данной статье показано как получить АЧХ и ФЧХ для любого фильтра с помощью Python имея передаточную функцию комплексной переменной s. Если вкратце, то для полчения АЧХ и ФЧХ нужно в выражение передаточной функции вместо s подставить j\omega и вычислить модуль и аргумент полученного комплексного выражения.

P.S. Должен заметить, что АЧХ и ФЧХ может изменятся в зависимости от параметров A и Q