【Scala】高速フーリエ変換を使って周期データをいじってみる #Scala

依存

// FFT用のライブラリ
libraryDependencies += "com.github.wendykierp" % "JTransforms" % "3.0"

// グラフ描画用のライブラリ
libraryDependencies += "com.github.sh0nk" %  "matplotlib4j" % "0.4.0"

正弦波をフーリエ変換してみる

import java.lang

import org.jtransforms.fft.DoubleFFT_1D
import com.github.sh0nk.matplotlib4j.Plot

object Main extends App {
  // １秒間のデータ数
  val sampleRate = 44100
  // 今回は882Hzの正弦波を使います
  val hertz = 882

  // 標本数(FFTで解析するには最大周波数の2倍以上の標本が必要)
  val bufferSize = 128

  val fft = new DoubleFFT_1D(bufferSize)

  // 空の配列を用意します
  val buffer: Array[Double] = Array.ofDim(bufferSize)
  // 882Hzの正弦波データ配列を生成します
  for (i <- 0 until bufferSize) {
    buffer(i) = periodSin(i, sampleRate, hertz)
  }

  // 正しく生成されたか確認してみましょう
  plot(buffer, "original", "time [s]", "amplitude", 1d / sampleRate)

  // フーリエ変換実行
  fft.realForward(buffer)

  // 変換後のbufferの中は: 奇数番 -> 実部, 偶数番目 -> 虚部　のようになっていて、それらの絶対値が対応する周波数の強度になっています
  val mag: Array[Double] = Array.ofDim(bufferSize / 2)
  for (i <- 0 until (bufferSize, 2)) {
    mag(i / 2) = Math.sqrt(buffer(i) * buffer(i) + buffer(i + 1) * buffer(i + 1))
  }

  // 横軸の周波数に注意して関係を確認してみましょう
  plot(mag, "strength", "frequency [Hz]", "strength", sampleRate / bufferSize)



  def periodSin(i: Int, sampleRate: Int, hertz: Int): Double = {
    Math.sin(Math.PI * 2 * i * hertz / sampleRate)
  }

  def plot(list: Array[Double], title: String = "title", xLabel: String, yLabel: String, dx: Double, dy: Double = 1d): Unit = {
    import java.util
    import collection.JavaConverters._
    val x: util.List[lang.Double] = list.indices.map(a => Double.box(a * dx)).asJava
    val y: util.List[lang.Double] = list.map(a => Double.box(a * dy)).toBuffer.asJava
    val plt = Plot.create
    plt.plot.add(x, y, ".")
    plt.xlabel(xLabel)
    plt.ylabel(yLabel)
    plt.title(title)
    plt.show()
  }

}

プロット図

変換前波形
変換後の周波数と強度の関係
拡大してみると800Hz前後がピークになっている

ポイント

標本数(FFTで解析するには最大周波数の2倍以上の標本が必要)
今回は882Hzの正弦波のため(44100 / 882 * 2 = 100)以上の標本があればいいのでbufferSize = 128とした
変換後のデータは連続した二つのデータがセットになっています(奇数番 -> そのセットの実部, 偶数番目 -> そのセットの虚部)
それぞれのセットは、そのセットの先頭から数えた数 * sampleRate / bufferSizeHzについてのデータにあたります。

正弦波を逆フーリエ変換してみる

  // 上のコードの続き

  fft.realInverse(buffer, true)
  // 逆変換後の結果
  plot(buffer, "after", "time [s]", "amplitude", 1d / sampleRate)

プロット図

元に戻っていることがわかる

周波数でフィルターしてみる

import java.lang

import org.jtransforms.fft.DoubleFFT_1D
import com.github.sh0nk.matplotlib4j.Plot

object Main extends App {
  // １秒間のデータ数
  val sampleRate = 44100
  // 今回は882Hzの正弦波を使います
  val hertz = 882
  val highCut = 2400

  // 標本数(FFTで解析するには最大周波数の2倍以上の標本が必要)
  val bufferSize = 512

  val fft = new DoubleFFT_1D(bufferSize)

  // 空の配列を用意します
  val buffer: Array[Double] = Array.ofDim(bufferSize)
  // 882Hzの正弦波データ配列を生成します
  for (i <- 0 until bufferSize) {
    // 2400Hzの正弦波を混ぜてみます
    buffer(i) = periodSin(i, sampleRate, hertz) + periodSin(i, sampleRate, highCut)
  }

  // 正しく生成されたか確認してみましょう
  plot(buffer, "original", "time [s]", "amplitude", 1d / sampleRate)

  // フーリエ変換実行
  fft.realForward(buffer)

  // 変換後のbufferの中は: 奇数番 -> 実部, 偶数番目 -> 虚部　のようになっていて、それらの絶対値が対応する周波数の強度になっています
  val mag: Array[Double] = Array.ofDim(bufferSize / 2)
  for (i <- 0 until (bufferSize, 2)) {
    mag(i / 2) = Math.sqrt(buffer(i) * buffer(i) + buffer(i + 1) * buffer(i + 1))
  }

  // 横軸の周波数に注意して関係を確認してみましょう
  plot(mag, "strength", "frequency [Hz]", "strength", sampleRate / bufferSize)



  // フーリエ変換後のbufferは連続した二つがセットになっていることに注意する
  // 2400Hzのデータは 2400 / (sampleRate / bufferSize) セットに相当する
  // 今回は大雑把に200Hz以上の音を全てカットする
  for (i <- 2000 / (sampleRate / bufferSize) until bufferSize) {
    buffer(i) = 0
  }


  fft.realInverse(buffer, true)
  // 逆変換後の結果
  plot(buffer, "after", "time [s]", "amplitude", 1d / sampleRate)


  def periodSin(i: Int, sampleRate: Int, hertz: Int): Double = {
    Math.sin(Math.PI * 2 * i * hertz / sampleRate)
  }

  def plot(list: Array[Double], title: String = "title", xLabel: String, yLabel: String, dx: Double, dy: Double = 1d): Unit = {
    import java.util
    import collection.JavaConverters._
    val x: util.List[lang.Double] = list.indices.map(a => Double.box(a * dx)).asJava
    val y: util.List[lang.Double] = list.map(a => Double.box(a * dy)).toBuffer.asJava
    val plt = Plot.create
    plt.plot.add(x, y, ".")
    plt.xlabel(xLabel)
    plt.ylabel(yLabel)
    plt.title(title)
    plt.show()
  }

プロット図

変換前波形
変換後の周波数と強度の関係
2000Hz以上をカットしてフーリエ逆変換

かなり雑にカットしたので完全に元の形にはなりませんがうまいことカットできたようです。
※緑がカット前、青がカット後

ポイント

変換後のデータ配列の消したい周波数に対応しているセットの値を0にした後、逆変換することで、その周波数をカットできる

終わりに

フーリエ変換はなんとなくで勉強したけど、初めて面白さがわかった気がします笑
今度は実際の音声データ(.wave)からノイズをカットしてみようかなと思ってます！