PinkTSP信号をJuliaで作る

Julia言語でPinkTSP信号を作るプログラムを書いた。とは言っても、以前書いたMatlabのものをそのまま移植しただけ。そもそも参考にしたのは以下の文献。

• 藤本卓也. “低域バンドでのSN比改善を目的としたTSP信号に関する検討.” 日本音響学会研究発表会講演論文集, pp.433–432. 1999年9月.
• 守谷直也, 金田豊. “Logarithmic TSP 信号を用いた高調波歪の検討.” 日本音響学会研究発表会講演論文集, pp.637–638. 2004年3月.

最低でも2秒間の信号（4回の同期加算）を44100Hz/16bitで欲しいときには

julia> pinktsp_generate(2.0, 44100, 16, 4);

とすると、pinktsp_2972msec_44100Hz_16bits_4times.wavというファイルができる。1回あたりの信号長は2秒よりも長い2の累乗の長さになるようにしてあるので、結果として2¹⁷÷44100=2.972秒となる（2の累乗の信号長にしておくとFFTが効率よく使えるので）。掃引方向を上向きにしたり（上記守谷文献参照）、スピーカ再生時にクリックノイズが出ないようにしたり、実務上の理由によりオリジナルの藤本文献とは違うことをいくつかしている。

この信号を使って録音したものをインパルス応答に戻すプログラムはまだ書いてない。また後日、時間があるときに。

コードは以下。

using SampledSignals, LibSndFile, Primes

function pinktsp_generate(duration, samprate=48000, nbits=32, reps=4)
    # duration = minimum length of the signal (seconds)
    # fs = number of samples per second (Hz)
    # nbits = number of bits per sample (bit)
    # reps = number of repeated measurements (times)

    N = nextpow2(ceil(Int, duration * samprate));
    m = prevprime(N >> 2) + 1;

    a = 2 * m * pi / ( (N/2) * log(N/2) );

    H = complex(zeros(N));
    H[1] = 1;
    H[2:convert(Int, N/2+1)] = exp.(1im * a .* (1:N/2) .* log.(1:N/2)) ./ sqrt.(1:N/2);
    H[convert(Int, N/2+2):N] = conj(H[convert(Int, N/2):-1:2]);
    h = real(ifft(H));

    # avoid click at the beginning and the end of the signal
    mi = indmin(abs.(h));
    h = [h[mi:end] ; h[1:mi-1]];
    h = reverse(h);

    # repeat
    hh = [repeat(h, outer=reps) ; zeros(length(h))];

    # write to an audio file
    hh = hh / maximum(abs.(hh)) * sqrt(1/2);   # peak at -3 dBFS
    if nbits==16
        hh2 = map(PCM16Sample, hh);
    #elseif nbits==24
    #    hh2 = map(PCM24Sample, hh);
    elseif nbits==32
        hh2 = map(PCM32Sample, hh);
    else
        nbits = 32;
        hh2 = map(PCM32Sample, hh);
        error("unsupported bitrate (falls back to 32bits)");
    end
    buf = SampleBuf(hh2, samprate);
    filename = @sprintf("pinktsp_%dmsec_%dHz_%dbits_%dtimes.wav", round(N/samprate*1000), samprate, nbits, reps);
    save(filename, buf);
end

実行環境は

julia> versioninfo()
Julia Version 0.6.2
Commit d386e40c17 (2017-12-13 18:08 UTC)
Platform Info:
  OS: macOS (x86_64-apple-darwin14.5.0)
  CPU: Intel(R) Core(TM) i7-7920HQ CPU @ 3.10GHz
  WORD_SIZE: 64
  BLAS: libopenblas (USE64BITINT DYNAMIC_ARCH NO_AFFINITY Prescott)
  LAPACK: libopenblas64_
  LIBM: libopenlibm
  LLVM: libLLVM-3.9.1 (ORCJIT, broadwell)