-- amstutz-1978-symmetric.lua -- -- Pierre Amstutz: -- Elliptic Approximation and Elliptic Filter Design on Small Computers, -- IEEE Transactions on Circuits and Systems (Volume: 25, Issue 12, December 1978). -- -- 2017-09-29 typed into Lua, H. Henkel pi, sin, cos = math.pi, math.sin, math.cos exp, log, sqrt = math.exp, math.log, math.sqrt abs, int = math.abs, math.floor fmt = string.format --********************************************************************** print("*:*:* SYMMETRICAL ELLIPTIC FILTERS *:*:*") --********************************************************************** b, c, d, e, f = {}, {}, {}, {}, {} dbn = log(10) / 10 -- dbn = 0.23025851 -- pi = 3.1415927 -- 2000*---------------------------------------------------------------- repeat repeat repeat repeat print("STOPBAND EDGE (KHZ) ?") fs = io.read("*n") print("PASSBAND EDGE (KHZ) ?") fp = io.read("*n") until abs(fs - fp) > 0 print("NUMBER OF PEAKS (1-15) ?") n = io.read("*n") until n > 0 m = 2 * n + 1 -- fc = sqrt(fs * fp) r = fc + fc for k = 1, 2 do s = fs + fp for j = 1, 6 do p = sqrt(s * r) s = (s + r) / 2 if 1.0e8 * (s - p) < s then break end r = p end if k >= 2 then break end q = m / s r = abs(fs - fp) end q = q * s s = exp(-pi / q) y = s print(fmt("CRITICAL Q %-16.9g", q / (4 * (1 - s) * s^n))) -- 2240* print("STOPBAND REJECTION (DB) ?") s = io.read("*n") until s > 0 s = exp(s * dbn / 2) r = exp(pi * q) p = (log(1 + (s * s - 1) / (r / 4 + 1 / r)^2)) / dbn print(fmt("PASSBAND RIPPLE (DB) %-16.9g", p)) r = r / (2 * (s + sqrt(s * s - 1))) r = log(r + sqrt(r * r + 1)) / (2 * q) r = sin(r) / cos(r) w = r print(fmt("3 DB ABOUT (KHZ) %-16.9g", fp + (fs - fp) / (1 + fc / (fp * r * r)))) print("NOMINAL RESISTANCE (OHM) ?") r = io.read("*n") until r > 0 -- 2380* z = y e[n] = w w = w * w for j = 1, m - 1 do f[j] = 1 end k = 1 for j = 1, 1024 do f[k] = f[k] * (1 - z) / (1 + z) if k >= m - 1 then z = z * y x = ((1 - z) / (1 + z))^2 e[n] = e[n] * (w + x) / (1 + w * x) k = 0 end z = z * y if z < 0.25e-18 then break end k = k + 1 end for j = 1, n do f[j] = f[j] * f[m - j] f[m - j] = f[j] end -- 3000*---------------------------------------------------------------- for j = 1, n do d[j] = f[2 * j] * (1 - f[j]^4) / f[j] b[j] = e[n] * f[j] end c[1] = 1 / b[n] for j = 1, n - 1 do c[j + 1] = (c[j] - b[n - j]) / (1 + c[j] * b[n - j]) e[n - j] = e[n + 1 - j] + e[n] * d[j] / (1 + b[j] * b[j]) end -- 4000*---------------------------------------------------------------- for j = 1, n do b[j] = ((1 + c[j] * c[j]) * e[j] / d[j] - c[j] / f[j]) / 2 c[j] = c[j] * f[j] d[j] = f[j] * f[j] end b[n + 1] = b[n] c[n + 1] = c[n] d[n + 1] = d[n] -- 5000*---------------------------------------------------------------- if n ~= 1 then for l = 1, 2 do for k = l + 2, n + 1, 2 do for j = l, k - 2, 2 do y = c[j] - c[k] z = 1 / (y / (b[j] * (d[k] - d[j])) - 1) b[k] = (b[k] - b[j]) * z * z - b[j] * (1 + z + z) c[k] = y * z end end end -- 6000*---------------------------------------------------------------- s = b[n] / b[n + 1] - 1 end q = 0.0005 / (pi * fc) p = q * r q = q / r if fs >= fp then -- 6050* print(" *:* LOW-PASS FILTER *:*") for j = 1, n do c[j] = q * c[j] d[j] = q * b[j] * d[j] b[j] = p / b[j] f[j] = fc / f[j] end c[n + 1] = q * c[n + 1] else -- 6140* print(" *:* HIGH-PASS FILTER *:*") for j = 1, n do c[j] = q / c[j] d[j] = q / (b[j] * d[j]) b[j] = p * b[j] f[j] = fc * f[j] end c[n + 1] = q / c[n + 1] end -- 6220* print(" KILOHERTZ FARAD HENRY") for j = 1, n, 2 do print(fmt(" %-17.8g", c[j])) print(fmt(" %-15.8g %2d %-17.8g%-17.8g", f[j], j, d[j], b[j])) end print(fmt(" %-17.8g", c[n + 1])) if n == 1 then exit(0) end l = (int((n + 1) / 2)) * 2 k = m - 1 - l for j = l + 2, m - 1, 2 do print(fmt(" %-15.8g %2d %-17.8g%-17.8g", f[k], k, d[k], b[k])) print(fmt(" %-17.8g", c[k])) k = k - 2 end print(fmt("PRECISION TEST %-16.9g", s))