インパルス応答測定用信号の生成

Angelo Farina氏のExponential Sine Sweep(ESS)をJuliaで実装しています http://pcfarina.eng.unipr.it/Public/Presentations/AES122-Farina.pdf

実行すると

  • ESS.wav(測定用信号)
  • InverseFilter.wav(逆フィルタ)
  • IR.wav(ESSと逆フィルタを畳み込んで作られたインパルス)
    の3つが作られ、それらのプロットが表示されます。
In [8]:
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)
Out[8]:
17.723884952769538
In [9]:
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)
Out[9]:

SampleBuf display requires javascript

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.

In [3]:
# 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")
Out[3]:
PyObject Text(0.5,1,'Amplitude envelope for inverse filter')
In [4]:
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)
Out[4]:

SampleBuf display requires javascript

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.

In [5]:
IR = conv(ESS, InverseFilter)
IR /= maximum(abs.(IR))
plot(IR)
grid()
title("Impulse response")
Buffer_IR = SampleBuf(IR, FS)
Out[5]:

SampleBuf display requires javascript

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.

In [6]:
save("ESS.wav", Buffer_ESS)
save("InverseFilter.wav", Buffer_InverseFilter)
save("IR.wav", Buffer_IR)