忘備録-備忘録

技術的な備忘録

広告

※このエリアは、60日間投稿が無い場合に表示されます。記事を投稿すると、表示されなくなります。

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.*;
  10. float gscale = 10; // 縦軸の倍率
  11. float goffset = 30; //オフセット
  12. float fmag = 2; //横軸の倍率
  13. int samplingrate = 1000; //サンプリングレート
  14. int displayxsize = 1000; //ディスプレイの横幅
  15. int displayysize = 600; //ディスプレイの高さ
  16. String comport = "COM5"; //シリアルポート
  17. int comspeed = 115200; //通信速度
  18. //fft data size 2^N
  19. int DATASIZE = 1024;
  20. Serial myPort;
  21. float val;
  22. float [] x = new float[DATASIZE]; //グラフ用データ x:時間
  23. float [] y1 = new float[DATASIZE]; //グラフ用データ y:電圧
  24. float [] y2 = new float[DATASIZE]; //グラフ用データ y:電圧
  25. //FFT用
  26. int fftdatasize = DATASIZE;
  27. FFT fft = new FFT(fftdatasize);
  28. float[] fftdata1r = new float[fftdatasize]; //data1実数部
  29. float[] fftdata1i = new float[fftdatasize]; //data1虚数部
  30. float[] fftdata2r = new float[fftdatasize]; //data2実数部
  31. float[] fftdata2i = new float[fftdatasize]; //data2虚数部
  32. float[] fftdata1 = new float[fftdatasize]; //data1絶対値
  33. float[] fftdata2 = new float[fftdatasize]; //data2絶対値
  34. void setup()
  35. {
  36.   frameRate(50); //20msec毎にグラフを更新
  37.   size(displayxsize, displayysize);
  38.   myPort = new Serial(this, comport, comspeed);
  39.   textSize(20);
  40.   //グラフ用データの初期化
  41.   for (int i = 0; i < x.length; i++) {
  42.     x[i] = (float)i * (float)(displayxsize-100)/(float)DATASIZE*fmag;
  43.     y1[i] = 0;
  44.     y2[i] = 0;
  45.   }
  46. }
  47. //縦軸の数値の加工 対数化
  48. float dispfunc(float data)
  49. {
  50.   float ans;
  51.   ans = gscale *( log(data) + goffset);
  52.   return ans;
  53. }
  54. void draw()
  55. {
  56.   background(0);
  57.   //グラフエリアの描画
  58.   fill(255);
  59.   text((int)(samplingrate/fmag) + "Hz", displayxsize - 80, displayysize-10);
  60.   float t = DATASIZE/1000;
  61.   text("0Hz", 30, displayysize-10);
  62.   fill(0, 255, 0);
  63.   for (int i = 0; i < fftdata1r.length; i++)
  64.   {
  65.     float[] wnd = fft.getWindow();
  66.     fftdata1[i] = (fftdata1r[i]*fftdata1r[i]+fftdata1i[i]*fftdata1i[i]);
  67.     fftdata2[i] = (fftdata2r[i]*fftdata2r[i]+fftdata2i[i]*fftdata2i[i]);
  68.     fftdata1r[i] = y1[i] * wnd[i];
  69.     fftdata2r[i] = y2[i] * wnd[i];
  70.     fftdata1i[i] = fftdata2i[i]=0;
  71.   }
  72.   fft.fft(fftdata1r, fftdata1i);
  73.   fft.fft(fftdata2r, fftdata2i);
  74.   //グラフのプロット
  75.   pushMatrix();
  76.   translate(50, displayysize-50);
  77.   scale(1, -1);
  78.   fill(0);
  79.   stroke(255);
  80.   rect(0, 0, displayxsize-100, displayysize-94);
  81.   stroke(0, 255, 0);
  82.   strokeWeight(1);
  83.   for (int i = 0; i < x.length/fmag - 1; i++) {
  84.     line(x[i], dispfunc(fftdata1[i]), x[i+1], dispfunc(fftdata1[i+1]));
  85.   }
  86.   stroke(255, 0, 0);
  87.   for (int i = 0; i < x.length/fmag - 1; i++) {
  88.     line(x[i], dispfunc(fftdata2[i]), x[i+1], dispfunc(fftdata2[i+1]));
  89.   }
  90.   popMatrix();
  91. }
  92. int bufsize = 512;
  93. char[] buf = new char[bufsize];
  94. int i = 0;
  95. void serialEvent(Serial myPort) {
  96.   char c = char(myPort.read());
  97.   if ((i<bufsize)&&(c!='\n'))
  98.   {
  99.     buf[i] = c;
  100.     i++;
  101.   } else
  102.   {
  103.     buf[i]='\0';
  104.     i=0;
  105.     String s = new String(buf);
  106.     float data[] = float(s.split(",", 0));
  107.     //グラフ用データの更新
  108.     for (int i = 0; i < y1.length - 1; i++) {
  109.       y1[i] = y1[i+1];
  110.       y2[i] = y2[i+1];
  111.     }
  112.     y1[y1.length - 1] = data[0];
  113.     y2[y2.length - 1] = data[1];
  114.   }
  115. }

プログラム 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.  */
  22. public class FFT {
  23.   int n, m;
  24.   // Lookup tables. Only need to recompute when size of FFT changes.
  25.   float[] cos;
  26.   float[] sin;
  27.   float[] window;
  28.   public FFT(int n) {
  29.     this.n = n;
  30.     this.m = (int)(log(n) / log(2));
  31.     // Make sure n is a power of 2
  32.     if (n != (1<<m))
  33.       throw new RuntimeException("FFT length must be power of 2");
  34.     // precompute tables
  35.     cos = new float[n/2];
  36.     sin = new float[n/2];
  37.     for (int i=0; i<n/2; i++) {
  38.       cos[i] = cos(-2*PI*i/n);
  39.       sin[i] = sin(-2*PI*i/n);
  40.     }
  41.     makeWindow();
  42.   }
  43.   protected void makeWindow() {
  44.     // Make a blackman window:
  45.     // w(n)=0.42-0.5cos{(2*PI*n)/(N-1)}+0.08cos{(4*PI*n)/(N-1)};
  46.     window = new float[n];
  47.     for (int i = 0; i < window.length; i++)
  48.       window[i] = 0.42 - 0.5 * cos(2*PI*i/(n-1))
  49.         + 0.08 * cos(4*PI*i/(n-1));
  50.   }
  51.   public float[] getWindow() {
  52.     return window;
  53.   }
  54.   /***************************************************************
  55.    * fft.c
  56.    * Douglas L. Jones
  57.    * University of Illinois at Urbana-Champaign
  58.    * January 19, 1992
  59.    * http://cnx.rice.edu/content/m12016/latest/
  60.    *
  61.    * fft: in-place radix-2 DIT DFT of a complex input
  62.    *
  63.    * input:
  64.    * n: length of FFT: must be a power of two
  65.    * m: n = 2**m
  66.    * input/output
  67.    * x: float array of length n with real part of data
  68.    * y: float array of length n with imag part of data
  69.    *
  70.    * Permission to copy and use this program is granted
  71.    * as long as this header is included.
  72.    ****************************************************************/
  73.   public void fft(float[] x, float[] y)
  74.   {
  75.     int i, j, k, n1, n2, a;
  76.     float c, s, e, t1, t2;
  77.     // Bit-reverse
  78.     j = 0;
  79.     n2 = n/2;
  80.     for (i=1; i < n - 1; i++) {
  81.       n1 = n2;
  82.       while ( j >= n1 ) {
  83.         j = j - n1;
  84.         n1 = n1/2;
  85.       }
  86.       j = j + n1;
  87.       if (i < j) {
  88.         t1 = x[i];
  89.         x[i] = x[j];
  90.         x[j] = t1;
  91.         t1 = y[i];
  92.         y[i] = y[j];
  93.         y[j] = t1;
  94.       }
  95.     }
  96.     // FFT
  97.     n1 = 0;
  98.     n2 = 1;
  99.     for (i=0; i < m; i++) {
  100.       n1 = n2;
  101.       n2 = n2 + n2;
  102.       a = 0;
  103.       for (j=0; j < n1; j++) {
  104.         c = cos[a];
  105.         s = sin[a];
  106.         a += 1 << (m-i-1);
  107.         for (k=j; k < n; k=k+n2) {
  108.           t1 = c*x[k+n1] - s*y[k+n1];
  109.           t2 = s*x[k+n1] + c*y[k+n1];
  110.           x[k+n1] = x[k] - t1;
  111.           y[k+n1] = y[k] - t2;
  112.           x[k] = x[k] + t1;
  113.           y[k] = y[k] + t2;
  114.         }
  115.       }
  116.     }
  117.   }
  118. }
コメント

シリアルプロッタを作る

2016-01-27 18:02:26 | processing
Arduinoに追加されたシリアルプロッタは手軽に使えて動作の確認も簡単です。
自動的に、縦軸のレンジが自動で変化するのはお手軽に使えて良いのですが、サンプリングレートが高くなると表示が変化しすぎて目が回るような感じになります。
縦軸、横軸のレンジを決めて動作するプログラムがほしくなりました。手軽に使えるProcessing言語を利用して作ってみました。ひな型になるプログラムが「Arduino シリアル通信 その1:processing」にあったのでそのプログラムを手直ししたものです。

シリアル受信データ
2つの数字をカンマ区切りで一定間隔に出力したものです。シリアルポートから受信します。
-0.369629,0.001532
-0.450195,0.002479
-0.546875,0.003281
-0.650879,0.004101
-0.796875,0.004864

動作画面


プログラムリスト
  1. // 参考サイト http://shirotsuku.sakura.ne.jp/blog/?p=192
  2. import processing.serial.*;
  3. int displaytime = 4; //表示時間
  4. float gscale = 1; // 縦軸の倍率
  5. int samplingrate = 1000; //サンプリングレート
  6. int displayxsize = 1200; //ディスプレイの横幅
  7. int displayysize = 800; //ディスプレイの高さ
  8. String comport = "COM5"; //シリアルポート
  9. int comspeed = 115200; //通信速度
  10. int DATASIZE = samplingrate * displaytime;
  11. //int lf = 0x0d; // ASCII linefeed
  12. //シリアルポートを使用する
  13. Serial myPort;
  14. float val;
  15. float [] x = new float[DATASIZE]; //グラフ用データ x:時間
  16. float [] y1 = new float[DATASIZE]; //グラフ用データ y:電圧
  17. float [] y2 = new float[DATASIZE]; //グラフ用データ y:電圧
  18. void setup()
  19. {
  20.   frameRate(50); //20msec毎にグラフを更新
  21.   size(displayxsize, displayysize);
  22.   myPort = new Serial(this, comport, comspeed);
  23.   //myPort.bufferUntil(lf);
  24.   textSize(20);
  25.   //グラフ用データの初期化
  26.   for (int i = 0; i < x.length; i++) {
  27.     x[i] = (float)i * (float)(displayxsize-100)/(float)DATASIZE;
  28.     y1[i] = 0;
  29.     y2[i] = 0;
  30.   }
  31. }
  32. void draw()
  33. {
  34.   background(0);
  35.   //グラフエリアの描画
  36.   fill(255);
  37.   // text("-1.6V",10,310);
  38.   // text("+1.6V",10,54);
  39.   text("0sec", displayxsize - 70, displayysize-10);
  40.   float t = DATASIZE/1000;
  41.   text("-" + t + "sec", 30, displayysize-10);
  42.   fill(0, 255, 0);
  43.   //text(val/255.0*5,560,100);
  44.   //グラフのプロット
  45.   pushMatrix();
  46.   translate(50, displayysize-50);
  47.   scale(1, -1);
  48.   fill(0);
  49.   stroke(255);
  50.   rect(0, 0, displayxsize-100, displayysize-94);
  51.   stroke(0, 255, 0);
  52.   strokeWeight(1);
  53.   for (int i = 0; i < x.length - 1; i++) {
  54.     line(x[i], y1[i], x[i+1], y1[i+1]);
  55.   }
  56.   stroke(255, 0, 0);
  57.   for (int i = 0; i < x.length - 1; i++) {
  58.     line(x[i], y2[i], x[i+1], y2[i+1]);
  59.   }
  60.   popMatrix();
  61. }
  62. int bufsize = 512;
  63. char[] buf = new char[bufsize];
  64. int i = 0;
  65. void serialEvent(Serial myPort) {
  66.   char c = char(myPort.read());
  67.   //if(c=='\r') return;
  68.   if ((i<bufsize)&&(c!='\n'))
  69.   {
  70.     buf[i] = c;
  71.     i++;
  72.   } else
  73.   {
  74.     buf[i]=0;
  75.     i=0;
  76.     String s = new String(buf);
  77.     float data[] = float(s.split(",", 0));
  78.     //グラフ用データの更新
  79.     for (int i = 0; i < y1.length - 1; i++) {
  80.       y1[i] = y1[i+1];
  81.       y2[i] = y2[i+1];
  82.     }
  83.     // 0〜(displayysize-94)の間で変化するように
  84.     y1[y1.length - 1] = (data[0]+0)*(displayysize-94)*gscale/2 + (displayysize-94)/2;
  85.     // 0〜(displayysize-94)の間で変化するように
  86.     y2[y2.length - 1] = (data[1]+0)*(displayysize-94)*gscale/2 + (displayysize-94)/2;
  87.     //print(data[0]);
  88.     //print(" ");
  89.     //println(s);
  90.   }
  91. }
コメント