1
2
3
4
5
6
7 import numpy as np
8 import matplotlib.pyplot as plt
9
10 def pmplot(omegas, frequencieslist, newfrequencieslist, numberofgrids, \
11 errorlist, iterationcount,amplitude,h):
12 normalizedfrequencies=omegas/np.pi
13 subplotheight=4
14 subplotwidth=10
15 plt.figure(figsize=(subplotwidth,((subplotheight+1)*iterationcount)))
16 for k in range (iterationcount):
17 freqaxis = np.array(frequencieslist[k])/(1.0*numberofgrids)
18 newfreqaxis = np.array(newfrequencieslist[k])/(1.0*numberofgrids)
19 errorarray=np.array(errorlist[k])
20 frequenciesarray=np.array(frequencieslist[k])
21 freqerror=np.take(errorarray, frequenciesarray)
22 newfrequenciesarray=np.array(newfrequencieslist[k])
23 newfreqerror=np.take(errorarray, newfrequenciesarray)
24 plt.subplot2grid((subplotheight*iterationcount, 1), \
25 (subplotheight*k, 0), rowspan=subplotheight)
26 plt.plot(normalizedfrequencies, errorlist[k], label='error')
27 plt.plot(freqaxis, freqerror, 'ks', markerfacecolor='none', \
28 markeredgecolor='red', label='previous extrema')
29 plt.plot(newfreqaxis, newfreqerror, 'o', markerfacecolor='none', \
30 markeredgecolor='blue', label='new extrema')
31 plt.xlabel(r'$\omega$/$\pi$')
32 plt.ylabel('Error')
33 plt.legend(loc='center left', bbox_to_anchor=(1, 0.5))
34 plt.axis(xmin=0, xmax=1)
35 plt.gca().set_title("Iteration {} error".format(k+1), size=10)
36 plt.tight_layout()
37 plt.savefig('pm.png')
38 plt.show()
39 plt.plot(normalizedfrequencies,amplitude)
40 plt.axis(xmin=0, xmax=1)
41 plt.xlabel(r'$\omega$/$\pi$')
42 plt.ylabel(r'$A(\omega)$')
43 plt.gca().set_title(r'A($\omega$)',size=10)
44 plt.savefig('pma.png')
45 plt.show()
46 plt.stem(h)
47 plt.xlabel('n')
48 plt.ylabel('h(n)')
49 plt.gca().set_title('h(n)',size=10)
50 plt.savefig('pmh.png')
51 plt.show()
52
53
54 def pmalgorithm(taplength,passbandfrequency,stopbandfrequency,passbandweight, \
55 stopbandweight,griddensity,smallnumber,maxiterationcount):
56 pi = np.pi
57 flag = 0
58 numberofextrema = ((taplength-1)//2)+2
59 frequencieslist=[]
60 newfrequencieslist=[]
61 errorlist=[]
62 numberofgrids = griddensity*numberofextrema
63 middlefrequency = (passbandfrequency+stopbandfrequency)/2
64 gridabscissa = range(numberofgrids+1)
65 gridabscissa = np.array(gridabscissa)
66 omegas = np.array(np.multiply(gridabscissa, pi/numberofgrids));
67 weights = np.zeros(numberofgrids+1);
68 for i in range(0, numberofgrids+1):
69 if omegas[i] <= passbandfrequency:
70 weights[i] = passbandweight
71 elif omegas[i] >= stopbandfrequency:
72 weights[i] = stopbandweight
73 desiredresponse = np.zeros(numberofgrids+1)
74 for i in range(0, numberofgrids+1):
75 if omegas[i] <= middlefrequency:
76 desiredresponse[i] = 1
77 nontransitionbandfrequency = []
78 for i in range(0, numberofgrids+1):
79 if weights[i] > 0:
80 nontransitionbandfrequency.append(i)
81 if i!=numberofgrids:
82 if weights[i] > 0 and weights[i+1] == 0:
83 bandedge1 = i
84 if i != numberofgrids:
85 if weights[i] == 0 and weights[i+1] != 0:
86 bandedge2 = i+1
87 nontransitionbandfrequencies = np.array(nontransitionbandfrequency)
88 nontransitionbandfrequencieslength = len(nontransitionbandfrequencies)
89 lastnontransitionbandfrequencyindex = nontransitionbandfrequencieslength-1
90 lastfrequency = nontransitionbandfrequencies[lastnontransitionbandfrequencyindex]
91 frequenciesindexes = np.linspace(0, lastnontransitionbandfrequencyindex, \
92 numberofextrema)
93 frequenciesindexes = np.round(frequenciesindexes, 0)
94 frequenciesindexes = frequenciesindexes.astype(int)
95 frequencies = np.take(nontransitionbandfrequencies, frequenciesindexes)
96 frequencieslastindex = numberofextrema-1
97 signarray = np.empty(numberofextrema, int)
98 for i in range(0, numberofextrema):
99 signarray[i] = (-1)**i
100 desiredresponseatextrema = np.array(numberofextrema)
101 iterationcount = 0
102 amplitude = np.zeros(len(gridabscissa))
103 relativedifference = np.zeros(maxiterationcount)
104 newfrequencies0 = np.empty(numberofextrema, dtype=np.int32)
105 while(1):
106 iterationcount = iterationcount+1
107 frequencieswithoutlast = np.delete(frequencies, frequencieslastindex)
108 cosinesoffrequencies = np.cos(frequencies*np.pi/(numberofgrids))
109 cosinesoffrequencieswithoutlast = np.cos(frequencieswithoutlast*pi/numberofgrids)
110 desiredresponseatextrema = np.take(desiredresponse, frequencies)
111 weightsatextrema = np.take(weights, frequencies)
112 cosineomegas = np.array(len(omegas))
113 cosineomegas = np.cos(omegas)
114 lagrangeweights = np.zeros(numberofextrema)
115 lagrangeweights2 = np.zeros(numberofextrema-1)
116 for i in range(0, numberofextrema):
117 lagrangeweights[i]=1
118 for j in range(0, numberofextrema):
119 if j != i:
120 lagrangeweights[i] = lagrangeweights[i]*(cosinesoffrequencies[i] - \
121 cosinesoffrequencies[j])
122 lagrangeweights[i] = 1/lagrangeweights[i]
123 gamma0 = lagrangeweights
124 deviationnumerator = sum(np.divide(np.multiply(gamma0, desiredresponseatextrema), weightsatextrema))
125 deviationdenominator = sum(np.multiply(gamma0, np.divide(signarray, weightsatextrema)))
126 deviation = deviationnumerator/deviationdenominator
127 rextrema = np.subtract(desiredresponseatextrema, (np.divide(np.multiply(signarray, deviation), \
128 weightsatextrema)))
129 rm1extrema=np.delete(rextrema, len(rextrema)-1)
130 for i in range(0, numberofextrema-1):
131 lagrangeweights2[i]=1
132 for j in range(0,numberofextrema-1):
133 if j != i:
134 lagrangeweights2[i] = lagrangeweights2[i]*(cosinesoffrequencieswithoutlast[i] - \
135 cosinesoffrequencieswithoutlast[j])
136 lagrangeweights2[i] = 1/lagrangeweights2[i]
137 beta0 = lagrangeweights2
138 for i in range(0, len(gridabscissa)):
139 if not(np.isin(gridabscissa[i], frequencies)):
140 amplitudenumeratorsum = 0.0
141 amplitudedenominatorsum = 0.0
142 for j in range(0,len(frequencies)-1):
143 amplitudenumeratorsum = amplitudenumeratorsum + \
144 (rm1extrema[j]*((beta0[j])/(cosineomegas[i]-cosinesoffrequencies[j])))
145 amplitudedenominatorsum = amplitudedenominatorsum + \
146 (beta0[j]/(cosineomegas[i]-cosinesoffrequencies[j]))
147 amplitude[i] = amplitudenumeratorsum/amplitudedenominatorsum
148 else:
149 for j in range(0, len(frequencies)):
150 if (frequencies[j] == i):
151 amplitude[i] = rextrema[j]
152 err = np.multiply(np.subtract(amplitude, desiredresponse), weights)
153 j = 0
154 for i in range(1, len(err)-1):
155 if abs(err[i]) > abs(err[i-1]) and abs(err[i]) > abs(err[i+1]):
156 newfrequencies0[j] = i;
157 j = j+1
158 if bandedge1 in newfrequencies0:
159 pass
160 else:
161 newfrequencies0[j] = bandedge1
162 j = j+1
163 if bandedge2 in newfrequencies0:
164 pass
165 else:
166 newfrequencies0[j] = bandedge2
167 j = j+1
168 if abs(err[0]) > abs(err[len(err)-1]):
169 newfrequencies0[j] = 0
170 j = j+1
171 if j < numberofextrema:
172 newfrequencies0[j] = len(err)-1
173 else:
174 newfrequencies0[j] = len(err)-1
175 j = j+1
176 if j < numberofextrema:
177 newfrequencies0[j] = 0
178 newfrequencies = np.sort(newfrequencies0)
179 relativedifference[iterationcount] = (max(abs(err))-abs(deviation))/abs(deviation)
180 frequencieslist.append([])
181 newfrequencieslist.append([])
182 errorlist.append([])
183 freqlist=frequencies.tolist()
184 newfreqlist=newfrequencies.tolist()
185 errlist=err.tolist()
186 for j in range (len(frequencies)):
187 frequencieslist[iterationcount-1].append(freqlist[j])
188 newfrequencieslist[iterationcount-1].append(newfreqlist[j])
189 for j in range (len(err)):
190 errorlist[iterationcount-1].append(errlist[j])
191 if relativedifference[iterationcount] < smallnumber:
192 print("I have converged at iteration: ", iterationcount)
193 flag=1
194 break
195 frequencies=newfrequencies
196 if iterationcount == maxiterationcount:
197 break
198 M = taplength//2
199 mlinspace = np.linspace(0, M, M+1)
200 m =np.array (mlinspace)
201 frequenciesomegas = np.take(omegas, frequencies)
202 momegas=np.outer(frequenciesomegas, m)
203 cosmomegas=np.cos(momegas)
204 frequenciesweights = np.take(weights, frequencies)
205 lastcolumn = np.divide(signarray, frequenciesweights)
206 lastcolumn = lastcolumn.reshape(len(frequencies), 1)
207 matrix0 = np.hstack((cosmomegas, lastcolumn))
208 desired = np.take(desiredresponse, frequencies)
209 desiredtranspose = np.transpose(desired)
210 inversematrix = np.linalg.inv(matrix0)
211 x = np.dot(inversematrix, desiredtranspose)
212 a = x[0:len(x)-1]
213 first = a[-1:0:-1]/2
214 second = np.array(a[0])
215 third = a[1:M+2]/2
216 h = np.append(first, second)
217 h = np.append(h, third)
218 if flag == 1:
219 pmplot(omegas, frequencieslist, newfrequencieslist, numberofgrids, \
220 errorlist, iterationcount, amplitude, h)
221
222
223 def pmtest():
224 griddensity = 20
225 smallnumber=1e-7
226 passbandweight = 1
227 stopbandweight = 6
228 maxiterationcount=20
229 passbandfrequency = 0.32*np.pi
230 stopbandfrequency = 0.40*np.pi
231 taplength = 39
232 pmalgorithm(taplength, passbandfrequency, stopbandfrequency, passbandweight, \
233 stopbandweight, griddensity, smallnumber, maxiterationcount)
234
235
236 pmtest()