Fixed to work with Python3

[cweickhmann]

The following code will run in Python 3.

# ELLIPTIC FILTER DESIGN FOR A CLASS OF GENERALIZED
#               HALFBAND FILTERS
#
#               BY RASHID ANSARI
#
#                INTEPRETED BY
#                  R. J. TONG
# 
# CITATION:
#    R. Ansari, "Elliptic filter design for a class of generalized halfband filters," 
#    in IEEE Transactions on Acoustics, Speech, and Signal Processing, vol. 33, no. 5, pp. 1146-1150, Oct 1985.
#    doi: 10.1109/TASSP.1985.1164709
#    URL: http://ieeexplore.ieee.org/stamp/stamp.jsp?tp=&arnumber=1164709&isnumber=26193
#
# DATE:
#    2/26/2017
# 
# INPUTS:
#    fp = pass band freq. 0<fp<0.4999 i.e 0.44
#    dp = pass band ripple i.e. 0.0001
#    ds = stop band ripple i.e. 0.1
#
# Modified for Python3+ by C. Weickhmann <c.weickhmann@gsi.de>

import matplotlib.pyplot as plt
import numpy as np

fp = 0.495
dp = 0.001
ds = 0.001

wp = fp * np.pi
ws = np.pi - wp

print(f"wp = {wp/np.pi:.6f}")
print(f"ws = {ws/np.pi:.6f}")

d1 = ds
d2 = np.sqrt(2.0*dp - np.power(dp,2.0))
d = min([d1, d2])/2.0
d_prime = 1.0 - np.sqrt(1.0 - np.power(d,2.0))

print(f"d1 = {d1:f}, d2 = {d2:f}, d = {d:f}, d' = {d_prime:e}")

k = np.power(np.tan(wp/2.0), 2.0)
k_prime = np.sqrt(1.0 - np.power(k, 2.0))

print(f"k = {k:e}, k' = {k_prime:e}")

delta = (1.0 - np.power(d, 2.0)) / np.power(d, 2.0)
p = 0.5 * (1.0 - np.sqrt(k_prime)) / (1.0 + np.sqrt(k_prime))
q = p + 2.0*np.power(p, 5.0) + 15.0*np.power(p, 9.0) + 150.0*np.power(p, 13.0)

print(f"p = {p:e}, q = {q:e}")

N = int(np.ceil(2.0*np.log(4.0*delta)/-np.log(q)))
if N % 2 == 0:
	print('N is even, adding +1')
	N += 1

L = int((N - 1)/2)

print(f"N = {N:d}, L = {L:.1f}, delta = {delta:e}")

omegas = np.zeros(L)
m_max = 10
for n in range(L):
	i = n + 1
	num_sigma = 0.0
	# print('Calculating numerator...')
	for m in range(0, m_max):
		A = (-1.0)**m
		B = np.power(q, m*(m+1.0))
		C = np.sin((2.0*m + 1.0)*np.pi*i/N)
		num_sigma += A*B*C
		# pad = ' '*m
		# print('%s->%e'%(pad,num_sigma))

	den_sigma = 0.0
	# print('Calculating denominator...')
	for m in range(1, m_max):
		A = (-1.0)**m
		B = np.power(q, m*m)
		C = np.cos(2.0*m*np.pi*i/N)
		den_sigma += A*B*C
		# pad = ' '*(m-1)
		# print('%s->%e'%(pad,den_sigma))

	omega_i = 2.0*np.power(q, 0.25) * num_sigma / (1.0 + 2.0*den_sigma)
	omegas[n] = omega_i
	print(f"omega[{n:d}] = {omega_i:f}")

omegas_2 = omegas * omegas
ris = np.sqrt((1.0 - k*omegas_2)*(1.0 - omegas_2/k))
cos_thetas = (-1.0)**((np.arange(L)+1.0)+1.0) * ris/(1.0 + omegas**2.0)
ais = (1.0-cos_thetas)/(1.0+cos_thetas)
print(f"ai = {str(ais):s}")

h0_ai = ais[ais >= 1.0]
h1_ai = ais[ais < 1.0]

print(f"H0(z) = {str(h0_ai):s}")
print(f"H1(z) = {str(h1_ai):s}")

pts = 1000
F = np.arange(pts)/float(pts) * 0.5
z = np.exp(2.0j*np.pi*F)
z_2 = np.power(z, 2.0)

h0s = [(z_2 + ai) / (ai*z_2 + 1.0) for ai in h0_ai]
h1s = [(ai*z_2 + 1.0) / (z_2 + ai) for ai in h1_ai]

h0_total = np.ones(pts, dtype=np.complex64)
for h0 in h0s:
	h0_total *= h0

h1_total = np.ones(pts, dtype=np.complex64)
for h1 in h1s:
	h1_total *= h1

hs = 1.0/z*h0_total + h1_total

plt.plot(F, 20.0*np.log10(np.abs(hs)))
plt.show()