忘備録-備忘録

技術的な備忘録

FFTアナライザを作る

2016-01-28 12:38:02 | processing
「シリアルプロッタを作る」でプロッタを作成したので、フーリエ変換を行うプログラムも作ってみました。使用している言語はProcessingです。
FFT変換プログラムはFFT.javaのものを使わせてもらいました。

動作画面


プログラム RS232CgraphFFT03.pde
  1. // シリアルFFTアナライザ
  2. // シリアル入力の値をフーリエ変換する
  3. //
  4. // 参考サイト http://shirotsuku.sakura.ne.jp/blog/?p=192
  5. // http://hp.vector.co.jp/authors/VA046927/fft4gjava.html
  6. // https://www.ee.columbia.edu/~ronw/code/MEAPsoft/doc/html/FFT_8java-source.html
  7. //
  8. //
  9. import processing.serial.*;
  11. float gscale = 10; // 縦軸の倍率
  12. float goffset = 30; //オフセット
  13. float fmag = 2; //横軸の倍率
  15. int samplingrate = 1000; //サンプリングレート
  17. int displayxsize = 1000; //ディスプレイの横幅
  18. int displayysize = 600; //ディスプレイの高さ
  20. String comport = "COM5"; //シリアルポート
  21. int comspeed = 115200; //通信速度
  23. //fft data size 2^N
  24. int DATASIZE = 1024;
  26. Serial myPort;
  27. float val;
  28. float [] x = new float[DATASIZE]; //グラフ用データ x：時間
  29. float [] y1 = new float[DATASIZE]; //グラフ用データ y：電圧
  30. float [] y2 = new float[DATASIZE]; //グラフ用データ y：電圧
  32. //FFT用
  33. int fftdatasize = DATASIZE;
  34. FFT fft = new FFT(fftdatasize);
  35. float[] fftdata1r = new float[fftdatasize]; //data1実数部
  36. float[] fftdata1i = new float[fftdatasize]; //data1虚数部
  37. float[] fftdata2r = new float[fftdatasize]; //data2実数部
  38. float[] fftdata2i = new float[fftdatasize]; //data2虚数部
  40. float[] fftdata1 = new float[fftdatasize]; //data1絶対値
  41. float[] fftdata2 = new float[fftdatasize]; //data2絶対値
  43. void setup()
  44. {
  45.   frameRate(50); //20msec毎にグラフを更新
  46.   size(displayxsize, displayysize);
  47.   myPort = new Serial(this, comport, comspeed);
  48.   textSize(20);
  50.   //グラフ用データの初期化
  51.   for (int i = 0; i < x.length; i++) {
  52.     x[i] = (float)i * (float)(displayxsize-100)/(float)DATASIZE*fmag;
  53.     y1[i] = 0;
  54.     y2[i] = 0;
  55.   }
  56. }
  58. //縦軸の数値の加工 対数化
  59. float dispfunc(float data)
  60. {
  61.   float ans;
  62.   ans = gscale *( log(data) + goffset);
  63.   return ans;
  64. }
  66. void draw()
  67. {
  68.   background(0);
  70.   //グラフエリアの描画
  71.   fill(255);
  72.   text((int)(samplingrate/fmag) + "Hz", displayxsize - 80, displayysize-10);
  73.   float t = DATASIZE/1000;
  74.   text("0Hz", 30, displayysize-10);
  75.   fill(0, 255, 0);
  76.   for (int i = 0; i < fftdata1r.length; i++)
  77.   {
  78.     float[] wnd = fft.getWindow();
  79.     fftdata1[i] = (fftdata1r[i]*fftdata1r[i]+fftdata1i[i]*fftdata1i[i]);
  80.     fftdata2[i] = (fftdata2r[i]*fftdata2r[i]+fftdata2i[i]*fftdata2i[i]);
  81.     fftdata1r[i] = y1[i] * wnd[i];
  82.     fftdata2r[i] = y2[i] * wnd[i];
  83.     fftdata1i[i] = fftdata2i[i]=0;
  84.   }
  85.   fft.fft(fftdata1r, fftdata1i);
  86.   fft.fft(fftdata2r, fftdata2i);
  87.   //グラフのプロット
  88.   pushMatrix();
  89.   translate(50, displayysize-50);
  90.   scale(1, -1);
  91.   fill(0);
  92.   stroke(255);
  93.   rect(0, 0, displayxsize-100, displayysize-94);
  94.   stroke(0, 255, 0);
  95.   strokeWeight(1);
  96.   for (int i = 0; i < x.length/fmag - 1; i++) {
  97.     line(x[i], dispfunc(fftdata1[i]), x[i+1], dispfunc(fftdata1[i+1]));
  98.   }
  99.   stroke(255, 0, 0);
  100.   for (int i = 0; i < x.length/fmag - 1; i++) {
  101.     line(x[i], dispfunc(fftdata2[i]), x[i+1], dispfunc(fftdata2[i+1]));
  102.   }
  103.   popMatrix();
  104. }
  106. int bufsize = 512;
  107. char[] buf = new char[bufsize];
  108. int i = 0;
  110. void serialEvent(Serial myPort) {
  111.   char c = char(myPort.read());
  113.   if ((i<bufsize)&&(c!='\n'))
  114.   {
  115.     buf[i] = c;
  116.     i++;
  117.   } else
  118.   {
  119.     buf[i]='\0';
  120.     i=0;
  121.     String s = new String(buf);
  122.     float data[] = float(s.split(",", 0));
  123.     //グラフ用データの更新
  124.     for (int i = 0; i < y1.length - 1; i++) {
  125.       y1[i] = y1[i+1];
  126.       y2[i] = y2[i+1];
  127.     }
  128.     y1[y1.length - 1] = data[0];
  129.     y2[y2.length - 1] = data[1];
  130.   }
  131. }
プログラム fft.pde
  1. /*
  2.  * Copyright 2006-2007 Columbia University.
  3.  *
  4.  * This file is part of MEAPsoft.
  5.  *
  6.  * MEAPsoft is free software; you can redistribute it and/or modify
  7.  * it under the terms of the GNU General Public License version 2 as
  8.  * published by the Free Software Foundation.
  9.  *
  10.  * MEAPsoft is distributed in the hope that it will be useful, but
  11.  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12.  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
  13.  * General Public License for more details.
  14.  *
  15.  * You should have received a copy of the GNU General Public License
  16.  * along with MEAPsoft; if not, write to the Free Software
  17.  * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA
  18.  * 02110-1301 USA
  19.  *
  20.  * See the file "COPYING" for the text of the license.
  21.  */
  23. public class FFT {
  25.   int n, m;
  27.   // Lookup tables. Only need to recompute when size of FFT changes.
  28.   float[] cos;
  29.   float[] sin;
  31.   float[] window;
  33.   public FFT(int n) {
  34.     this.n = n;
  35.     this.m = (int)(log(n) / log(2));
  37.     // Make sure n is a power of 2
  38.     if (n != (1<<m))
  39.       throw new RuntimeException("FFT length must be power of 2");
  41.     // precompute tables
  42.     cos = new float[n/2];
  43.     sin = new float[n/2];
  45.     for (int i=0; i<n/2; i++) {
  46.       cos[i] = cos(-2*PI*i/n);
  47.       sin[i] = sin(-2*PI*i/n);
  48.     }
  50.     makeWindow();
  51.   }
  53.   protected void makeWindow() {
  54.     // Make a blackman window:
  55.     // w(n)=0.42-0.5cos{(2*PI*n)/(N-1)}+0.08cos{(4*PI*n)/(N-1)};
  56.     window = new float[n];
  57.     for (int i = 0; i < window.length; i++)
  58.       window[i] = 0.42 - 0.5 * cos(2*PI*i/(n-1))
  59.         + 0.08 * cos(4*PI*i/(n-1));
  60.   }
  62.   public float[] getWindow() {
  63.     return window;
  64.   }
  67.   /***************************************************************
  68.    * fft.c
  69.    * Douglas L. Jones
  70.    * University of Illinois at Urbana-Champaign
  71.    * January 19, 1992
  72.    * http://cnx.rice.edu/content/m12016/latest/
  73.    *
  74.    * fft: in-place radix-2 DIT DFT of a complex input
  75.    *
  76.    * input:
  77.    * n: length of FFT: must be a power of two
  78.    * m: n = 2**m
  79.    * input/output
  80.    * x: float array of length n with real part of data
  81.    * y: float array of length n with imag part of data
  82.    *
  83.    * Permission to copy and use this program is granted
  84.    * as long as this header is included.
  85.    ****************************************************************/
  86.   public void fft(float[] x, float[] y)
  87.   {
  88.     int i, j, k, n1, n2, a;
  89.     float c, s, e, t1, t2;
  92.     // Bit-reverse
  93.     j = 0;
  94.     n2 = n/2;
  95.     for (i=1; i < n - 1; i++) {
  96.       n1 = n2;
  97.       while ( j >= n1 ) {
  98.         j = j - n1;
  99.         n1 = n1/2;
  100.       }
  101.       j = j + n1;
  103.       if (i < j) {
  104.         t1 = x[i];
  105.         x[i] = x[j];
  106.         x[j] = t1;
  107.         t1 = y[i];
  108.         y[i] = y[j];
  109.         y[j] = t1;
  110.       }
  111.     }
  113.     // FFT
  114.     n1 = 0;
  115.     n2 = 1;
  117.     for (i=0; i < m; i++) {
  118.       n1 = n2;
  119.       n2 = n2 + n2;
  120.       a = 0;
  122.       for (j=0; j < n1; j++) {
  123.         c = cos[a];
  124.         s = sin[a];
  125.         a += 1 << (m-i-1);
  127.         for (k=j; k < n; k=k+n2) {
  128.           t1 = c*x[k+n1] - s*y[k+n1];
  129.           t2 = s*x[k+n1] + c*y[k+n1];
  130.           x[k+n1] = x[k] - t1;
  131.           y[k+n1] = y[k] - t2;
  132.           x[k] = x[k] + t1;
  133.           y[k] = y[k] + t2;
  134.         }
  135.       }
  136.     }
  137.   }
  138. }


最新の画像［もっと見る］

コメントを投稿