Angelo Farina氏のExponential Sine Sweep(ESS)をJuliaで実装しています http://pcfarina.eng.unipr.it/Public/Presentations/AES122-Farina.pdf
実行すると
の3つが作られ、それらのプロットが表示されます。
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)
To enable for the whole notebook select "Trust Notebook" from the "File" menu. You can also trust this cell by re-running it. You may also need to re-run `using SampledSignals` if the module is not yet loaded in the Julia kernel, or `SampledSignals.embed_javascript()` if the Julia module is loaded but the javascript isn't initialized.
# 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")
PyObject Text(0.5,1,'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)
To enable for the whole notebook select "Trust Notebook" from the "File" menu. You can also trust this cell by re-running it. You may also need to re-run `using SampledSignals` if the module is not yet loaded in the Julia kernel, or `SampledSignals.embed_javascript()` if the Julia module is loaded but the javascript isn't initialized.
IR = conv(ESS, InverseFilter)
IR /= maximum(abs.(IR))
plot(IR)
grid()
title("Impulse response")
Buffer_IR = SampleBuf(IR, FS)
To enable for the whole notebook select "Trust Notebook" from the "File" menu. You can also trust this cell by re-running it. You may also need to re-run `using SampledSignals` if the module is not yet loaded in the Julia kernel, or `SampledSignals.embed_javascript()` if the Julia module is loaded but the javascript isn't initialized.
save("ESS.wav", Buffer_ESS)
save("InverseFilter.wav", Buffer_InverseFilter)
save("IR.wav", Buffer_IR)