In [1]:
get_ipython().ast_node_interactivity = 'all' import matplotlib.pyplot as plt import numpy as np import matplotlib from scipy.signal import lfilter from scipy.optimize import * matplotlib.rcParams['figure.dpi'] = 150 def fft(x): x = np.fft.fft(x) x = np.abs(x) x = x / np.linalg.norm(x) x = x[0:5000] return x
In [2]:
class EMA_Filter: def __init__(self, alpha): self.alpha = alpha self.prev = None def __call__(self, sample): if self.prev is None: self.prev = sample return sample res = self.alpha * sample + (1 - self.alpha) * self.prev self.prev = res return res
In [3]:
impulse = np.zeros(10_000) impulse[5000] = 1 _ = plt.plot(fft(impulse), label="Flat") for alpha in [0.8, 0.6, 0.5]: f = EMA_Filter(alpha) _ = plt.plot(fft(list(map(f, impulse))), label=str(alpha)) _ = plt.legend()
Out:
Let’s try to replicate EMA_Filter(0.6) with IIR.
In [4]:
f = EMA_Filter(0.6) target = fft(list(map(f, impulse))) def prep_ab(x): split_ratio = x[0] l = len(x) h = int(l * split_ratio) x = x[1:] a = [1, *x[:h]] b = x[h:] return a, b def fitness(params): a, b = prep_ab(params) f = fft(lfilter(b, a, impulse)) return np.sum((target - f) ** 2) n = 2 x = minimize(fitness, [0.5] + [0.1] * n, bounds=[(0, 1)] + [(-1, 1)] * n) round(x.fun, 5) [round(x, 2) for x in x.x] a, b = prep_ab(x.x) print(f"alpha = {a}") print(f"beta = {b}") _ = plt.plot(fft(lfilter(b, a, impulse)), label="replicated") _ = plt.plot(target, label="target") _ = plt.legend()
Out [4]:
0.0
Out [4]:
[0.5, -0.4, 0.1]
Out:
alpha = [1, -0.39999887119697985] beta = [0.1]
Out:
In [8]:
def prep_ab(x): l = len(x) h = l // 2 a = [1, *x[:h]] b = x[h:] return a, b def safe_log10(x, eps=1e-10): result = np.where(x > eps, x, -10) np.log10(result, out=result, where=result > 0) return result def fitness(params): a, b = prep_ab(params) x = lfilter(b, a, impulse) x = np.clip(x, -1, 1) f = fft(x) if np.isnan(np.min(f)): return 9999 error = 0 error -= np.sum(f[:400] ** 2) error += np.sum(f[450:1000]) return error n = 15 #x = minimize(fitness, [0.1] * (n * 2), bounds=[(-2, 2)] * (n * 2)) #x = basinhopping(fitness, [0.1] * (n * 2)) x = dual_annealing(fitness, bounds=[(-1, 1)] * (n * 2)) #x = differential_evolution(fitness, bounds=[(-1, 1)] * (n * 2), maxiter=10) x.nfev x.message round(x.fun, 5) x = [round(x, 3) for x in x.x] a, b = prep_ab(x) print(f"alpha = {a}") print(f"beta = {b}") x = lfilter(b, a, impulse) x = np.clip(x, -1, 1) x = fft(x) plt.plot(x)
Out [8]:
62698
Out [8]:
['Maximum number of iteration reached']
Out [8]:
-0.56853
Out:
alpha = [1, -0.481, 0.24, -0.527, -0.28, -0.446, 0.123, 0.968, -0.489, 0.674, -0.924, -0.041, -0.781, 0.506, -0.458, -0.192] beta = [1.0, 0.945, 0.59, 0.746, 0.601, 0.836, -0.201, 0.163, -0.476, 0.353, -0.663, 0.832, -0.243, -0.076, 0.332]
Out [8]:
[]
Out:
In [6]:
x = 20 * safe_log10(x / np.max(x)) plt.plot(x)
Out [6]:
[]
Out:
In [7]:
a = [1, 0.3, 0.25, -0.1] b = [0.05, 0.01] plt.plot(fft(lfilter(b, a, impulse)))
Out [7]:
[]
Out:
