using LibSndFile
using SampledSignals
using PyPlot
rc("figure",figsize =([16,4]), ) #pyplotのサイズ

FS = 48000;#サンプルレート
T = 1;#信号長（秒）
NUMSAMPLES = FS * T;
FreqBegin = 20.0;#測定終了周波数（Hz）
FreqEnd = 24000.0;#測定終了周波数（Hz）

ESS = zeros(NUMSAMPLES);#Exponential Sine Sweep（ESS）配列
InverseFilter = zeros(NUMSAMPLES);#逆フィルタ配列　
AmplitudeEnvelope = zeros(NUMSAMPLES);#逆フィルタ用の音量フェード

alpha = (2pi * FreqBegin * T) / log(FreqEnd / FreqBegin)

for i in 1:NUMSAMPLES
    t = (i - 1) / FS
    ESS[i] = sin(alpha * (exp((t / T) * log(FreqEnd / FreqBegin)) - 1.0));
end;

#リンギング除去の為にESS終端の直近ゼロ地点までトリミングする
for i in NUMSAMPLES:-1:1
    if(abs.(ESS[i]) < 0.05)# ≒0
        break;
    else
        ESS[i] = 0.0
    end
end;

plot(ESS)
grid()
title("Exponential sine sweep")
Buffer_ESS = SampleBuf(ESS, FS)

# for i in 1:NUMSAMPLES
#     AmplitudeEnvelope[i] = 10.0 ^ ((-6.0 * log2(FreqEnd / FreqBegin)) * ((i - 1) / NUMSAMPLES)) * 0.05
# end
for i in 1:NUMSAMPLES
    AmplitudeEnvelope[i] = 10.0 ^ (((-6.0 * log2(FreqEnd / FreqBegin)) * ((i - 1) / NUMSAMPLES)) * 0.05)
end
plot(AmplitudeEnvelope)
grid()
title("Amplitude envelope for inverse filter")

n = NUMSAMPLES;
for i in 1:NUMSAMPLES
    InverseFilter[i] = ESS[n] * AmplitudeEnvelope[i];
    n-=1;
end
plot(InverseFilter)
grid()
title("Inverse filter")
Buffer_InverseFilter = SampleBuf(InverseFilter, FS)

IR = conv(ESS, InverseFilter)
IR /= maximum(abs.(IR))
plot(IR)
grid()
title("Impulse response")
Buffer_IR = SampleBuf(IR, FS)

save("ESS.wav", Buffer_ESS)
save("InverseFilter.wav", Buffer_InverseFilter)
save("IR.wav", Buffer_IR)