import numpy as np
import matplotlib.pyplot as plt
import math
import wave

# デジタル信号 f[i] から DFT 係数の実数成分 A[k] と虚数成分 B[k] を計算する関数
# 演習で作成したソースからコピーする
def DFT(A,B,f):
?

# DFT係数の実数成分 A[k] と虚数成分 B[k] から絶対値 C_abs[k] と偏角 C_arg[k] を計算する関数
# 演習で作成したソースからコピーする
def GetAbsArg(C_abs,C_arg,A,B):
?

# Waveファイル読み込み関数
def read_wavedata(name):
    print('read: '+name)
    with wave.open(name,'r') as wavein:
        print(wavein.getparams())
        rate = wavein.getframerate()
        nframes = wavein.getnframes()
        data = wavein.readframes(wavein.getnframes())
    return np.frombuffer(data, dtype=np.int16),rate,nframes

x,fs,frames = read_wavedata('DFT_2_1.wav')

print('')
print(f'サンプリング周波数 fs = {fs} Hz' )

# DFT_2_1.wav の start_sec 秒時点から dulation 秒間だけ音声信号を切り出してディジタル信号 f[i] を作成する
# 上手くグラフが表示されない時は適宜 start_sec を調整する
start_sec = 0.5
dulation = 0.05

# 周期 (= 信号長)
N = int(dulation*fs)
print(f'周期 N = {N} 点' )

# 基本周波数
f1 = fs/N
print(f'基本周波数 f1 = {f1} Hz')

# 基本周波数を元にインデックス k をヘルツに対応させる
freq=[ f1*k for k in range(0,int(N/2)) ]

# 音声信号の切り出し
start_pos = int( start_sec*fs )
end_pos = start_pos + N
f = x[start_pos:end_pos]

A=np.zeros(N) # DFT係数 C[k] の実数成分
B=np.zeros(N) # C[k] の虚数成分
C_abs=np.zeros(N) # C[k] の絶対値
C_arg=np.zeros(N) # C[k] の偏角

DFT(A,B,f)
GetAbsArg(C_abs,C_arg,A,B)

print('')
print('振幅成分')
plt.plot(freq,C_abs[:int(N/2)])
plt.xlabel('Hz')
plt.show()

print('位相成分')
plt.plot(freq,C_arg[:int(N/2)])
plt.xlabel('Hz')
plt.show()