HowToDoIt

#!/usr/bin/env python3
# -*- coding: utf-8 -*- 
""" 
Created on Wed May 16 21:52:15 2018 
@author:  Viktor Signaievskyi 
@blog:    HowToDoIt.com.ua 
@e-mail:  [email protected] 
@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();