import numpy as np
import math
import cmath
import matplotlib.pyplot as plt
from scipy.fft import fft, ifft

figsize=10
dpi=100

N=8

# 元の信号
def f_original(i) :
    return 5 + 6*np.cos(1*2*np.pi/8*i+np.pi/4) +4*np.cos(2*2*np.pi/8*i-np.pi/4) +2*np.cos(3*2*np.pi/8*i+np.pi/2) -6*np.cos(4*2*np.pi/8*i)

# 実フーリエ級数展開
def fourier(i,DFT,N) :
    out = DFT[0]
    M = int(N/2)
    if N%2==0:
        for k in range(1,M):
            out += 2*abs(DFT[k])*math.cos(k*2*math.pi/N*i+cmath.phase(DFT[k]))
        out += DFT[M]*math.cos(M*2*math.pi/N*i)
    else:
        for k in range(1,M+1):
            out += 2*abs(DFT[k])*math.cos(k*2*math.pi/N*i+cmath.phase(DFT[k]))
    return out.real

def draw_f():
    plt.xticks(np.arange(0,2*N,N))
    plt.axvline(N,color="gray",linestyle='dashed')
    plt.axhline(0,color="gray",linestyle='dashed')
    plt.xlabel(r'$i$')
    plt.ylabel(r'$f[i]$',rotation=0,labelpad=0,loc="top")
    plt.xlim(0,2*N-1)
    plt.ylim(-10,20)
    plt.show()
    
f = [f_original( i ) for i in range(N)]

DFT = fft(f)/N  # N で割っているのはテキストのDFTの定義とスケールを合わせるため

for k in range(N):
    print(f'{abs(DFT[k]):7.4f}, {cmath.phase(DFT[k])/np.pi:7.4f} : {DFT[k]} ')

f_restored = [fourier(i,DFT,N) for i in range(2*N)]

# 元の信号
plt.plot(range(2*N), np.r_[f,f],marker='o',color='blue')
plt.plot(range(2*N), np.r_[f,f],marker='o',markersize=5,linestyle='',color='black')
draw_f()

# 復元した信号
plt.plot(range(2*N), f_restored,marker='o',color='red')
plt.plot(range(2*N), f_restored,marker='o',markersize=5,linestyle='',color='black')
draw_f()