-
Notifications
You must be signed in to change notification settings - Fork 1
/
bpf_fir.py
executable file
·146 lines (117 loc) · 3.15 KB
/
bpf_fir.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
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
#!/usr/bin/python
import numpy as np
import scipy as sp
import scipy.signal as sg
import matplotlib.pyplot as plt
import pylab as pl
#first trying with given values
f_num=102
del_1=0.1
del_2=0.1
f_sample=90000
if f_num > 70:
m=f_num-70
else:
m=f_num
if (m*0.1-int(m*0.1)) > 0:
q_m=int(0.1*m)
else:
q_m=(0.1*m - 1)
r_m=m-10*q_m
B_l_m=2+0.7*q_m+2*r_m
B_h_m=B_l_m+5
# analog values for bandpass
omega_s1=(B_l_m-1)*1000.0
omega_p1=(B_l_m)*1000.0
omega_p2=(B_h_m)*1000.0
omega_s2=(B_h_m+1)*1000.0
#analog f_range array
f_analog_a=np.array([omega_s1,omega_p1,omega_p2,omega_s2],dtype='f')
#normalized digital frequencies
f_digital_a=(f_analog_a/f_sample)*2*np.pi
#Kaiser window parameters
del_omega1=f_digital_a[3]-f_digital_a[2]
del_omega2=f_digital_a[1]-f_digital_a[0]
del_omega=min(abs(del_omega1),abs(del_omega2))
A=-20*np.log10(del_1)
#Order
N_1=(A-8)/(2*2.285*del_omega)
N=np.ceil(N_1)
if(A<21):
alpha=0
elif(A<=50):
alpha=0.5842*(A-21)**0.4+0.07886(A-21)
else:
alpha=0.1102(A-8.7)
omega_c1=(f_digital_a[1]+f_digital_a[0])*0.5
omega_c2=(f_digital_a[3]+f_digital_a[2])*0.5
#Ideal BPF h[n]
iterable=((np.sin(omega_c2*k)-np.sin(omega_c1*k))/(np.pi*k) for k in range(int(-N),int(N+1)))
h_ideal=np.fromiter(iterable,float)
h_ideal[N]=((omega_c2-omega_c1)/np.pi)
beta=alpha/N
#Generating Kaiser window
h_kaiser=sg.kaiser(2*N+1,beta)
h_org=h_ideal*h_kaiser
print "Filter Coefficients:\n",h_org
#print "Hideal",h_ideal
#print m,q_m,r_m,B_l_m,B_h_m
#print "Analog frequencies",f_analog_a
#print "Digital frequencies",f_digital_a
#print "Equivalent Digital frequencies",f_eqv_analog_a
#print "Equivalent Analog lpf freq",f_eqv_analog_lpf_a
#print np.sqrt(D_1),np.sqrt(D_2),B
#print "order=",N
#print A_k
#print "Poles",poles_a
#print "Numerator",numer
#print "Denominator:",denom
#plotting poles
#plt.figure(1)
#plt.grid(True)
#plt.scatter(poles_a.real,poles_a.imag,marker='x')
#plt.legend()
nyq_rate=f_sample/2
#------------------------------------------------
# Plot the FIR filter coefficients.
#------------------------------------------------
plt.figure(1)
plt.plot(h_org, 'bo-', linewidth=2)
plt.title('Filter Coefficients (%d taps)' % (2*N+1))
plt.grid(True)
#Plot Frequency response
plt.figure(2)
plt.clf()
plt.grid(True)
w,h= sg.freqz(h_org)
plt.plot((w/np.pi)*nyq_rate, np.absolute(h), linewidth=2)
plt.xlabel('Frequency (Hz)')
plt.ylabel('Gain')
plt.title('Frequency Response')
plt.ylim(-0.05, 1.2)
# Upper inset plot.
ax1 = plt.axes([0.42, 0.6, .45, .25])
plt.plot((w/np.pi)*nyq_rate, np.absolute(h), linewidth=2)
plt.xlim(8000.0,13000.0)
plt.ylim(0.9, 1.2)
plt.grid(True)
# Lower inset plot
ax2 = plt.axes([0.42, 0.25, .45, .25])
plt.plot((w/np.pi)*nyq_rate, np.absolute(h), linewidth=2)
plt.xlim(13000.0, 20000.0)
plt.ylim(0.0, 0.11)
plt.grid(True)
plt.figure(3)
plt.grid(True)
h_Phase = pl.unwrap(np.arctan2(np.imag(h),np.real(h)))
plt.plot(w/max(w),h_Phase)
plt.ylabel('Phase (radians)')
plt.xlabel(r'Normalized Frequency (x$\pi$rad/sample)')
plt.title(r'Phase response')
#Stem Diagram
plt.figure(4)
y = pl.linspace(0,77,77)
plt.stem(y,h_org,linefmt='b-', markerfmt='bo', basefmt='r-')
plt.title('Filter Coefficients (%d taps)' % (2*N+1))
plt.grid(True)
plt.show()