4 スペクトル解析

4.1 パワースペクトルの計算例

（2019年1月28日 Julia 1.1に対応）

using FileIO: load
import LibSndFile
using FFTW
using Plots
pyplot();

# リコーダー吹鳴音の読み込み、冒頭から1.8秒分だけ使う
# （基音が約524Hzの吹奏音が入っている）
y = load("recorder.flac");
fs = y.samplerate;
y1 = y[1:Int(fs*1.8)];

# 2の累乗になるようにzero-padding
y2 = zeros(nextpow(2, length(y1)));   
y2[1:length(y1)] = y1;

# フーリエ変換して絶対値を計算し、2乗してパワーにする
Y2 = abs.(fft(y2)) / length(y2);
Y2 = Y2 .^ 2;

# 正の周波数成分と負の周波数成分をまとめる
# （DCからナイキストまでの範囲（正）だけを抽出し、
# 　捨ててしまった負の部分を補正するために2倍する）
Y3 = Y2[1 : Int(length(Y2)/2) + 1];
Y3[2:end-1] = Y3[2:end-1] * 2.0;

# デシベルに変換
Y4 = 10 * log10.(Y3);

# 各ビンの周波数を計算
f = range(0, stop=y.samplerate/2, length=length(Y3));

# ここまでで、fとY4に各ビンの周波数とパワーが入った

# グラフ描画
plot(f, Y4 .- maximum(Y4),   # 最大値が0 dBになるようにして描画
     xlim=(0, 15000), ylim=(-80, +3),   # X軸・Y軸の範囲指定
     grid=true,   # グリッドを表示
     legend=false,   # 凡例を表示しない
     title="Power Spectrum",   # タイトルや軸のラベルを付ける
     xlabel="Frequency (Hz)",
     ylabel="Power (dB)");

savefig("output.png");   # PNGファイルへの保存
savefig("output.pdf");   # PDFファイルへの保存

↓こんなのが出てくる。

図: リコーダー吹鳴音のパワースペクトル