フーリエ逆変換を用いたローパスフィルタ(Player10)

2008/09/1

こんにちは。
私は子供にせがまれてよく絵(落書き)を描きます。今日は、こいのぼりを描いてといわれて、描くと今度はこいのぼりがベビーカーに乗っているように描けといいます。さらにベルトを締めて、自分が押す姿を描けとも。
こいのぼりがベビーカーにベルトをしめて乗っかり、子供がそれを押している絵ができました。シュールです。

さて、前回フーリエ変換ができました。それで、本によると、それをちょろっと変えるとフーリエ逆変換ができるらしいです。
波形→フーリエ変換で、周波数成分に分解。
周波数成分→フーリエ逆変換で、波形にもどす。
ということができるようになるわけです。ここで、周波数成分に分解したときに、特定の周波数以上、または以下を消すことで、ローパス・ハイパス・バンドパスフィルターと似たようなことができます。
普通はそんなまわりくどいことをせずに、FIR・IIRフィルタというような遅延処理を使ってやるみたいです。

ただ、今回フーリエ逆変換ができて小躍りしてしまったので、とりあえずアップです。
→サンプルはこちら(要Player10)

遊び方なんですが、わかりやすいのが、周波数を高めに設定して、カットオフ周波数を低めにすることです。
ある一定の部分から急に音がこもったり変化したらOKです。
ノイズでやると、よりわかりやすいかもです。

プログラム的には、
フーリエ変換した後に出てくる実部と虚部のデータから、必要ない部分を0にします。
それをフーリエ逆変換して、実際のサンプル数でわると、波形データが出てきます。
この辺、本を見ながらやっただけなんで、間違ってるかもしれません。まあ、遊びなんで許してください、、。

FFTのソースです。

package
{
  import __AS3__.vec.Vector;

  public class FFT
  {
    public static const MODE_FFT:String = "mode_fft";
    public static const MODE_IFFT:String = "mode_ifft";

    //FFT
    public static function getFFT(ar:Vector.<Number>):Vector.<Number>
    {
      var n:uint = ar.length;
      //一時処理用
      var temp:Vector.<Number> = new Vector.<Number>(n);

      //虚数部を全て0にする
      var count:uint = n;
      while(count > 0)
      {
        temp[count - 1] = 0;
        count--;
      }

      //FFT。実数部と虚数部を得る
      var resultFFT:Vector.<Vector.<Number>> = calcFFT(ar, temp, MODE_FFT);

      //ベクトルの大きさを求める
      var i:uint;
      for(i = 0; i < n; i++)
      {
         temp[i] = Math.sqrt(Math.pow(resultFFT[0][i], 2) + Math.pow(resultFFT[1][i], 2));
      }

      return temp;
    }

    //IFFT。入力に実数部と虚数部が必要
    public static function getIFFT(ar:Vector.<Number>, ai:Vector.<Number>):Vector.<Number>
    {
      //IFFT。実数部と虚数部を得る
      var resultFFT:Vector.<Vector.<Number>> = calcFFT(ar, ai, MODE_IFFT);

      var i:uint;
      var n:uint = resultFFT.length;
      var real:Vector.<Number> = resultFFT[0];
      for(i = 0; i < n; i++)
      {
        real[i] = real[i] / n;
      }

      //実数部のみ返す
      return real;
    }

    //ローパスフィルタ。cutoffFrequencty: 0.0 ~ 0.5
    public static function lowPassFilter(ar:Vector.<Number>, cutoffFrequency:Number):Vector.<Number>
    {
      var n:uint = ar.length;
      //一時処理用
      var ai:Vector.<Number> = new Vector.<Number>(n);

      //虚数部を全て0にする
      var i:uint;
      for(i = 0; i < n; i++)
      {
        ai[i] = 0;
      }

      //FFT。実数部と虚数部を得る
      var resultFFT:Vector.<Vector.<Number>> = calcFFT(ar, ai, MODE_FFT);

      //不必要な部分をカット
      ar = resultFFT[0];
      ai = resultFFT[1];
      var min:uint;
      var max:uint;
      min = n * cutoffFrequency;
      max = n - min;
      for(i = 0; i < n; i++)
      {
        if(i > min && i < max)
        {
          ar[i] = 0;
          ai[i] = 0;
        }
      }

      //IFFT
      var output:Vector.<Number> = getIFFT(ar, ai);
      return output;
    }

    //arは実数部, aiは虚数部
    public static function calcFFT(ar:Vector.<Number>, ai:Vector.<Number>, mode:String):Vector.<Vector.<Number>>
    {
      var n:uint = ar.length; //データの数。必ず2の乗数であること
      var theta:Number = 2 * Math.PI / n;
      var m:int, mh:int, i:int, j:int, k:int, irev:int;
      var wr:Number, wi:Number, xr:Number, xi:Number;

      //scrambler
      i = 0;
      for (j = 1; j < n - 1; j++) {
        for (k = n >> 1; k > (i ^= k); k >>= 1);
        if (j < i) {
          xr = ar[j];
          xi = ai[j];
          ar[j] = ar[i];
          ai[j] = ai[i];
          ar[i] = xr;
          ai[i] = xi;
        }
      }
      for (mh = 1; (m = mh << 1) <= n; mh = m) {
        irev = 0;
        for (i = 0; i < n; i += m) {
          wr = Math.cos(theta * irev);
          wi = Math.sin(theta * irev);
          for (k = n >> 2; k > (irev ^= k); k >>= 1);
          for (j = i; j < mh + i; j++) {
            k = j + mh;
            xr = ar[j] - ar[k];
            xi = ai[j] - ai[k];
            ar[j] += ar[k];
            ai[j] += ai[k];
            if(mode == MODE_FFT)
            {
              ar[k] = wr * xr - wi * xi;
              ai[k] = wr * xi + wi * xr;
            }
            else if(mode == MODE_IFFT)
            {
              ar[k] = wr * xr + wi * xi;
              ai[k] = wr * xi - wi * xr;
            }

          }
        }
      }

      var output:Vector.<Vector.<Number>> = new Vector.<Vector.<Number>>(2);
      output[0] = ar;
      output[1] = ai;
      return output;
    }
  }
}

あと、またまた参考までに、クライアント側のソースです。

package
{
  import __AS3__.vec.Vector;

  import fl.controls.Button;
  import fl.controls.ComboBox;
  import fl.controls.Label;
  import fl.controls.Slider;
  import fl.events.SliderEvent;

  import flash.display.Graphics;
  import flash.display.Sprite;
  import flash.display.StageScaleMode;
  import flash.events.Event;
  import flash.events.MouseEvent;
  import flash.events.SampleDataEvent;
  import flash.media.Sound;
  import flash.media.SoundChannel;
  import flash.text.TextField;
  import flash.text.TextFieldAutoSize;
  import flash.text.TextFormat;

  [SWF(width="640", height="350", frameRate="30", backgroundColor="#b6c3be")]
  public class IFFTsample extends Sprite
  {
    private var graphBase:Sprite;
    private var graph:Sprite;
    private var asset:FFT_asset2;
    private var isPlay:Boolean;
    private var ch1:SoundChannel;
    private var sound:Sound;
    private var waveTypeCombo:ComboBox;
    public static const SINE:String = "sine";
    public static const SAW:String = "saw";
    public static const TRIANGLE:String = "triangle";
    public static const PULSE:String = "pulse";
    public static const NOISE:String = "noise";
    private var waveType:String = SINE;
    private var frequencySlider:Slider;
    private var frequency:uint = 440;
    private var cutoffFreqSlide:Slider;
    private var cutoffFrequency = 0.5;

    public function IFFTsample()
    {
      super();
      initialize();
    }

    private function initialize():void
    {
      stage.scaleMode = StageScaleMode.NO_SCALE;
      initializeGraph();
      initializeAsset();

      frequency = 440;
      isPlay = true;
      sound = new Sound();
      sound.addEventListener("sampleData", sampleDataComplete);
    }

    private function sampleDataComplete(e:SampleDataEvent):void
    {
      var input:Vector.<Number> = new Vector.<Number>();
      var output:Vector.<Number>;
      var sample:Number;
      //var frequency:uint = 440;
      var f0:Number = frequency / 44100;
      var f1:Number = 44100 / frequency;
      var PAI2:Number = 2 * Math.PI;
      for ( var c:int=0; c < 2048; c++ )
      {
        switch(waveType)
        {
          //サイン波
          case SINE:
          sample = Math.sin(PAI2 * f0 * Number(c + e.position));
          break;

          //矩形波
          case PULSE:
          sample = (Math.sin(PAI2 * f0 * Number(c + e.position)) < 0) ? 1 : -1;
          break;

          //のこぎり波
          case SAW:
          sample = 2 * f0 * ((c + e.position) % f1) - 1;
          break;

          //三角波
          case TRIANGLE:
          sample = (Math.sin(PAI2 * f0 * Number(c + e.position)) < 0) ?
              (4 * f0 * ((c + e.position) % (f1 / 2)) - 1) :
              (-4 * f0 * ((c + e.position) % (f1 / 2)) + 1);
          break;

          //ホワイトノイズ
          case NOISE:
          sample = Math.random() * 2 - 1;
          break;

          default: break;
        }
        input[c] = sample;
      }

      output = FFT.lowPassFilter(input, cutoffFrequency);
      for(c = 0; c < 2048; c++)
      {
        sample = output[c]
        if(sample > 1)
        {
          sample = 1;
        }
        else if(sample < -1)
        {
          sample = -1;
        }
        e.data.writeFloat(output[c]);
        e.data.writeFloat(output[c]);
      }
      output.splice(output.length / 4, output.length / 4 * 3);
      drawWaveForm(output, graph, 380, 280);
    }

    private function drawFFT(data:Vector.<Number>, target:Sprite, WIDTH:uint, HEIGHT:uint):void
    {
      var g:Graphics = target.graphics;
      var i:uint;
      var len:uint = data.length;
      var rate:Number = WIDTH / len;
      var step:Number = (rate < 1) ? (1 / rate) : 1;
      g.clear();
      g.lineStyle(1, 0xff0000, 1);
      g.moveTo(0, HEIGHT);
      for(i = 0; i < data.length; i += step)
      {
        g.moveTo(i * rate, HEIGHT);
        g.lineTo(i * rate, HEIGHT - HEIGHT * data[i] / 800);
      }
    }

    private function drawWaveForm(data:Vector.<Number>, target:Sprite, WIDTH:uint, HEIGHT:uint):void
    {
      var g:Graphics = target.graphics;
      var i:uint;
      var len:uint = data.length;
      var rate:Number = WIDTH / len;
      var step:Number = (rate < 1) ? (1 / rate) : 1;
      var HALF:uint = HEIGHT / 2;
      g.clear();
      g.lineStyle(1, 0xff0000, 1);
      g.moveTo(0, HALF);
      for(i = 0; i < data.length; i += step)
      {
        g.lineTo(i * rate, HALF - HALF * data[i] / 3600);
      }
    }

    private function initializeGraph():void
    {
      var g:Graphics;
      var devideNum:uint = 10;
      var i:uint;
      var devideWidth:Number;
      var HEIGHT:uint;
      var tf:TextField;
      var tft:TextFormat;
      var samplingRate:uint = 44.1;
      graphBase = new Sprite();
      addChild(graphBase);
      graph = new Sprite();
      graphBase.addChild(graph);
      g = graphBase.graphics;
      g.beginFill(0x000000, 1);
      g.drawRect(0,0,380,280);
      g.endFill();
      graphBase.x = 15;
      graphBase.y = 10;
    }

    private function initializeAsset():void
    {
      asset = new FFT_asset2();
      addChild(asset);

      var btn:Button = asset.startStop;
      btn.label = "スタート";
      btn.addEventListener(MouseEvent.MOUSE_DOWN, startStopHandler);

      initializeCombobox();
      initializeFrequencySlider();
      initializeCutoffFreqSlide();
    }

    private function initializeCombobox():void
    {
      waveTypeCombo = asset.combo;
      waveTypeCombo.addItem( { label: "サイン波", data:SINE } );
      waveTypeCombo.addItem( { label: "矩形波", data:PULSE } );
      waveTypeCombo.addItem( { label: "のこぎり波", data:SAW } );
      waveTypeCombo.addItem( { label: "三角波", data:TRIANGLE } );
      waveTypeCombo.addItem( { label: "ノイズ", data:NOISE } );
      waveTypeCombo.addEventListener(Event.CHANGE, waveTypeComboChangeHandler);
    }

    private function waveTypeComboChangeHandler(e:Event):void
    {
      waveType = waveTypeCombo.value;
    }

    private function initializeFrequencySlider():void
    {
      frequencySlider = asset.freqSlide;
      frequencySlider.width = 200;
      frequencySlider.tickInterval = 20;
      frequencySlider.maximum = 12000;
      frequencySlider.minimum = 20;
      frequencySlider.value = frequency;
      frequencySlider.addEventListener(SliderEvent.THUMB_DRAG, frequencySliderDrag);
      var label:Label = asset.freq_txt;
      label.text = frequency + "Hz";
    }

    private function frequencySliderDrag(e:SliderEvent):void
    {
      frequency = frequencySlider.value;
      var label:Label = asset.freq_txt;
      label.text = frequency + "Hz";
    }

    private function initializeCutoffFreqSlide():void
    {
      cutoffFreqSlide = asset.cutoffFreqSlide;
      cutoffFreqSlide.width = 200;
      cutoffFreqSlide.tickInterval = 1;
      cutoffFreqSlide.maximum = 1000;
      cutoffFreqSlide.minimum = 0;
      cutoffFreqSlide.value = cutoffFrequency * 2000;
      cutoffFreqSlide.addEventListener(SliderEvent.THUMB_DRAG, cutoffFreqSlideDrag);
      var label:Label = asset.cutoffFreq_txt;
      label.text = cutoffFrequency;
    }

    private function cutoffFreqSlideDrag(e:SliderEvent):void
    {
      cutoffFrequency = cutoffFreqSlide.value / 2000;
      var label:Label = asset.cutoffFreq_txt;
      label.text = cutoffFrequency;
    }

    private function startStopHandler(e:MouseEvent):void
    {
      var btn:Button = asset.startStop;
      if(isPlay)
      {
        btn.label = "ストップ";
      }
      else
      {
        btn.label = "スタート";
      }
      playStopSound();
      isPlay = !isPlay;
    }

    private function playStopSound():void
    {
      if(isPlay)
      {
        ch1 = sound.play();
      }
      else
      {
        ch1.stop();
      }
      var g:Graphics = graph.graphics;
      g.clear();
    }
  }
}

自作iPhoneアプリ 好評発売中!
フォルメモ - シンプルなフォルダつきメモ帳
ジッピー電卓 - 消費税や割引もサクサク計算!

ページトップへ戻る