#!/usr/bin/env python
# coding: utf-8

# # ITU-T G.227 複素数フィルタ → 実数フィルタ
# 
# まず、ITU-T G.227 https://www.itu.int/rec/T-REC-G.227-198811-I/en
# 
# 記載の伝達関数を周波数特性としてグラフ化。普通のフィルタのグラフと違ってyはゲインではなくロスなので注意
# 
# 
# $$ \frac{E}{2V} = 
# \frac{
#   18400 + 91238\ p^2 + 11638\ p^4 + p (67280 + 54050\ p^2)
# }{
#   400 + 4001\ p^2 + p^4 + p (36040 + 130\ p^2)
# }
# $$
# $$
# p = j \frac{f}{1000}
# $$

# In[1]:


import numpy as np
from matplotlib import pyplot as plt
from scipy import signal
import json

def g227(niq, size):
    f = np.linspace(0, niq, size)
    p = 1j * f / 1000
    numerator = 18400 + 91238*p**2 + 11638*p**4 + p*(67280 + 54050*p**2)
    denominator = 400 + 4001*p**2 + p**4 + p*(36040 + 130*p**2)
    loss = np.abs(numerator / denominator)
    dB = 20 * np.log10(loss)
    # adjust gain
    dB -= np.min(dB)
    loss = np.power(10, dB / 20)
    return f, loss


# In[2]:


sr = 44100
firlen = 513
xticks = [33, 50, 100, 500, 1000, 5000, 10000]

f_m, loss_m = g227(sr/2, 2**13)
f, loss = g227(sr/2, int(firlen/2))

print(np.min(20 * np.log10(loss)))
plt.figure(figsize=(8,10))
plt.title("Frequency response")
plt.xscale('log')
plt.grid(True)
plt.ylabel('Composite loss')
plt.xlabel('Frequency')
plt.ylim((0, 70))
#plt.xlim((33, 10000))
plt.xlim((1, sr/2))
plt.xticks(xticks, xticks)
plt.plot(f, 20 * np.log10(loss), label="length=%d" % int(firlen/2))
plt.plot(f_m, 20 * np.log10(loss_m), label="expected")
plt.legend(bbox_to_anchor=(1, 1), loc='upper right', borderaxespad=1, fontsize=11)

plt.show()


# # FIR
# 
# そのまま IFFT して有限インパルス応答 (FIR) を得る。フィルタにしたいのでロスの逆数をとってゲインにしている。

# In[3]:


fir = np.fft.fftshift(np.fft.ifft(1.0/loss, n=firlen))
print(len(fir))
x = np.linspace(-firlen/2, firlen/2, len(fir))
plt.figure(figsize=(15,5))
plt.title("Impulse response")
plt.grid(True)
plt.ylabel("Coeffs")
plt.xlabel("Samples")
#plt.xlim((-30,30))
plt.plot(x, np.real(fir), label="real")
plt.plot(x, np.imag(fir), label="imag")
plt.legend(bbox_to_anchor=(1, 1), loc='upper right', borderaxespad=1, fontsize=11)
plt.show()


# # FIR の周波数特性
# 
# 求めた FIR の周波数特性を評価する。freqz を使う。
# 

# In[4]:


w, h = signal.freqz(fir, worN=2**13)
x = w * sr * 1.0 / (2 * np.pi)
plt.figure(figsize=(10,10))
plt.title("Frequency response of impulse response (complex)")
plt.xscale('log')
plt.grid(True)
plt.ylabel('Gain')
plt.xlabel('Frequency')
plt.xlim((1, sr/2))
plt.ylim((-100, 10))
plt.xticks(xticks, xticks)

freqz_dB = 20 * np.log10(abs(h))
expected_dB = -20 * np.log10(loss_m)
plt.plot(x, freqz_dB, 'b', label="freqz")
plt.plot(f_m, expected_dB, 'r', label="expected")
plt.legend(bbox_to_anchor=(1, 1), loc='upper right', borderaxespad=1, fontsize=11)
print(len(h) / np.sum(np.abs(h)))

error = np.abs(freqz_dB- expected_dB)
print('Max Error', np.max(error), 'dB', 'at', x[np.argmax(error)], 'Hz')
print('MSE',  20*np.log10(np.sqrt(np.mean(np.power(10, error / 20)**2))), 'dB')

ax2 = plt.twinx()
angles = np.unwrap(np.angle(h))
ax2.plot(x, angles, 'g')
ax2.set_ylabel('Angle (radians)', color='g')

plt.show()


# # 実数フィルタ
# 
# 実数フィルタにする。
# 

# In[5]:


firr = np.copy(np.real(fir))
firr *= 2

w, h = signal.freqz(firr, worN=2**13)
x = w * sr * 1.0 / (2 * np.pi)

plt.figure(figsize=(10,10))
plt.title("Frequency response of impulse response (real)")
plt.xscale('log')
plt.grid(True)
plt.ylabel('Gain')
plt.xlabel('Frequency')
plt.xlim((1, sr/2))
plt.ylim((-100, 10))
plt.xticks(xticks, xticks)

freqz_dB = 20 * np.log10(abs(h))
expected_dB = -20 * np.log10(loss_m)
plt.plot(x, freqz_dB, 'b', label="freqz")
plt.plot(f_m, expected_dB, 'r', label="expected")
plt.legend(bbox_to_anchor=(1, 1), loc='upper right', borderaxespad=1, fontsize=11)
print(len(h) / np.sum(np.abs(h)))
error = np.abs(freqz_dB- expected_dB)
print('Max Error', np.max(error), 'dB', 'at', x[np.argmax(error)], 'Hz')
print('MSE',  20*np.log10(np.sqrt(np.mean(np.power(10, error / 20)**2))), 'dB')

ax2 = plt.twinx()
angles = np.unwrap(np.angle(h))
ax2.plot(x, angles, 'g')
ax2.set_ylabel('Angle (radians)', color='g')
plt.show()


# # JSON

# In[6]:


print(json.dumps(list(map(lambda x: (x.real, x.imag), fir))))


# In[7]:


print(json.dumps(list(map(lambda x: x.real, firr))))


# In[ ]:





# In[ ]: