-
Notifications
You must be signed in to change notification settings - Fork 0
/
iir2_freq.py
128 lines (106 loc) · 3.15 KB
/
iir2_freq.py
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
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
#
# get 2nd order iir filter frequency response
#
# frequency range is from FCL to FCH
# calcuation points are divided by BAND_NUM as log-scale
# output unit is dB. And normalized frequency by sampling frequency FS0
#
# Soem of following code is refer from id:aidiary's Blog page.
# Check version
# Python 2.7.12 on win32 (Windows version)
# numpy 1.9.2
# matplotlib 1.5.3
import numpy as np
from matplotlib import pyplot as plt
FS0=48000
FCL=50.0
FCH=1000.0
BAND_NUM=100
class IIR2_FREQ:
def __init__(self,fcl=FCL,fch=FCH,fs0=FS0,band_num=BAND_NUM,plt_pause0=False):
self.FCL=fcl
self.FCH=fch
self.FS0=fs0
self.BAND_NUM=band_num
self.plt_pause=plt_pause0
def H2(self,f, a, b):
nume = b[0] + b[1] * np.exp(-2j * np.pi * f) + b[2] * np.exp(-4j * np.pi * f)
deno = 1 + a[1] * np.exp(-2j * np.pi * f) + a[2] * np.exp(-4j * np.pi * f)
val = nume / deno
return np.sqrt(val.real ** 2 + val.imag ** 2)
def H0(self,a,b):
amp=[]
freq=[]
bands= np.zeros( self.BAND_NUM+1)
fcl=(self.FCL * 1.0) / self.FS0 # fc * 1.0 means convert fc to float
fch=(self.FCH * 1.0) / self.FS0 # fc * 1.0 means convert fc to float
delta1=np.power(fch/fcl, 1.0 / (self.BAND_NUM)) # Log Scale
bands[0]=fcl
#print "i,band = 0", bands[0] * FS0
for i in range(1, self.BAND_NUM+1):
bands[i]= bands[i-1] * delta1
#print "i,band =", i, bands[i] * self.FS0
for f in bands:
amp.append(self.H2(f, a, b))
#print f
return np.log10(amp) * 20 ,bands * self.FS0
def H0disp(self,a,b):
amp, freq= self.H0(a,b)
plt.clf()
plt.plot(freq,amp)
plt.xlabel('Hz')
plt.ylabel('dB')
if self.plt_pause :
plt.pause(1.0)
else:
plt.show()
def H0disp2(self,a1,b1,a2, b2):
amp1, freq1= self.H0(a1,b1)
amp2, freq2= self.H0(a2,b2)
plt.clf()
plt.plot(freq1,amp1, label='initial')
plt.plot(freq2,amp2, label='target')
plt.xlabel('Hz')
plt.ylabel('dB')
plt.legend()
if self.plt_pause :
plt.pause(5.0)
else:
plt.show()
def H0disp3(self,a0,b0,a1,b1,a2, b2):
amp0, freq0= self.H0(a0,b0)
amp1, freq1= self.H0(a1,b1)
amp2, freq2= self.H0(a2,b2)
plt.clf()
plt.plot(freq0,amp0, label='initial')
plt.plot(freq1,amp1, label='present')
plt.plot(freq2,amp2, label='target')
plt.xlabel('Hz')
plt.ylabel('dB')
plt.legend()
plt.title('frequency response')
if self.plt_pause :
plt.pause(5.0)
else:
plt.show()
def disp3(self,a0,b0,a1,b1,a2, b2):
print "+++ comparison of iir filter coefficient"
print " initial present target"
print "a[0]", a0[0],a1[0],a2[0]
print "a[1]", a0[1],a1[1],a2[1]
print "a[2]", a0[2],a1[2],a2[2]
print "b[0]", b0[0],b1[0],b2[0]
print "b[1]", b0[1],b1[1],b2[1]
print "b[2]", b0[2],b1[2],b2[2]
if __name__ == '__main__':
import iir1
# sample plot pattern 1
pi0= iir1.IIR1(fs0=48000)
a1, b1= pi0.CUT(180,-8,8)
pif0=IIR2_FREQ(fcl=20,fch=20000,fs0=48000,band_num=1024,plt_pause0=True)
pif0.H0disp(a1,b1)
# sample plot pattern 2
a1, b1= pi0.CUT(200,-6,8) # initial
a2, b2= pi0.CUT(180,-8,8) # target if known
pif1=IIR2_FREQ()
pif1.H0disp2(a1,b1,a2,b2)