統計ソフトRでFFTして遊ぼう -コンプレッサーの周波数特性-

数学関数をお手軽にFFTして遊べたら楽しいのになーとずいぶん前から思っていた.
なんでそんなこと考えてたかというと, 波形をコンプで潰したときに
周波数がどうなるのかなーていうのを数学的にコントロールしたかったから.
発想としては思いつくけど実装が大変そうでなかなか実行できずにいたけど,
統計ソフトRを用いると本当に気軽にそれが実行出来るということを知った.

かなりこちらの記事を参考にさせてもらいました.
Diaspar Journal: 1分で試すフーリエ変換（FFT）

まずR上で変数等の初期設定.

sampling = 4096            # FFTサンプルポイント数
n = 0:(sampling-1)         # ポイント数分のベクトル
samplefreq = 44100         # サンプリング周波数(Hz)
t = n/samplefreq           # 時間軸(s)
f = n*samplefreq/sampling  # 周波数軸(Hz)
wave = sin(100*2*pi*t)     # 変換したい波形 この場合は100Hzのサイン波

samplingやsamplefreq, waveは適当に変更可能.
周波数解像度はsamplefreq/samplingで, 今の場合だと10Hzくらい.

そして波形や周波数を描画. 簡単簡単.

plot(t, wave, type="l", xlim=c(0,0.02), ylim=c(-1,1))         # 時間領域でプロット
spec = abs(fft(wave))^2                       # フーリエ変換, 2乗してパワースペクトル
plot(f, spec, type="l", xlim = c(20,samplefreq/2), log="xy")  # 周波数領域でプロット

コンプレッサー関数はこんな感じ.
閾値(1/2)を超えたらレシオ(4:1)にしたがって線形に潰す(もちろんこれらも任意に変更可能).
波がマイナスの時は, absでプラスに持ってきて潰してからsignで反転.
もっといい方法があるかも.
ifelseでベクトルに対して条件分岐できるのがRの便利なところっぽい.

compline = function(x) (x+1.5)/4
comp = function(x) {
  ifelse(abs(x)<1/2, x, sign(x)*compline(abs(x)))
}

コンプしたものを描画. 奇数倍音が立ってるのがわかります.

plot(t, comp(wave), type="l", xlim=c(0,0.02), ylim=c(-1,1))           
compspec = abs(fft(comp(wave)))^2                               
plot(f, compspec, type="l", xlim = c(20,samplefreq/2), log="xy")