Irr Ema Visualization

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

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

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()

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()

0.0

[0.5, -0.4, 0.1]

alpha = [1, -0.39999887119697985]
beta  = [0.1]

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)

62698

['Maximum number of iteration reached']

-0.56853

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]

[<matplotlib.lines.Line2D at 0x7efd2b32d900>]

x = 20 * safe_log10(x / np.max(x))
plt.plot(x)

[<matplotlib.lines.Line2D at 0x7efd2b45d690>]

a = [1, 0.3, 0.25, -0.1]
b = [0.05, 0.01]

plt.plot(fft(lfilter(b, a, impulse)))

[<matplotlib.lines.Line2D at 0x7efd2b49fdf0>]