(高速フーリエ変換)(SSD1306)ふたたび、UNOに周波数分布を表示する。(入力ランダムウォーク)(FFT)(高速化)





x 過去ログをみよ

結果

![o_coq349.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/840352/94a85460-334c-6222-3c0f-7d0926cb8e08.jpeg)



プログラム



```c++



//
//attiny85-spectrum_32band.ino
//https://github.com/colonelwatch/attiny85-spectrum/
//
//Arduino_UNO_3_spectrum_32band_3
//
//Ver 1.0 いろいろ
//Ver 2.0 いろいろ
//Ver 3.0 高速化　表示をオリジナル　64の割り算をマクロ　転記を減らす
//


//インクルド
#include "Arduino.h"
#include <Wire.h>
#include "fix_fft.h"                  // 8-bit FFT library modified for Arduino


//定義(フーリエ変換)
char data_a[64];                                                 // used to store FFT input/output and past data
//unsigned long microseconds;                                    // used for timekeeping
int summ, avg;                                                   // used for DC bias elimination
//#define SAMPLING_FREQUENCY 15000  // Sampling frequency (Actual max measured frequency captured is half)
//const unsigned int sampling_period_us = round(1000000 * (2.0 / SAMPLING_FREQUENCY)); // Sampling period (doubled to account for overclock)


//定義(SSD1306)
#define MAX_PAGE                   (7)
#define MAX_COL                    (127)

#define COMMAND_MODE               0x80 // continuation bit is set!
#define DATA_MODE                  0x40

#define SET_COLUMN_ADDRESS         0x21 // takes two bytes, start address and end address of display data RAM
#define SET_PAGE_ADDRESS           0x22 // takes two bytes, start address and end address of display data RAM

#define SET_MEMORY_ADDRESSING_MODE 0x20 // takes one byte as given above
#define HORIZONTAL_ADDRESSING_MODE 0x00

#define SET_SEGMENT_REMAP_0        0xA0 // column address 0 is mapped to SEG0 (Reset)
#define SET_SEGMENT_REMAP_127      0xA1 // column address 127 is mapped to SEG0

#define SET_COMMON_REMAP_0         0xC0 // row address  0 is mapped to COM0 (Reset)
#define SET_COMMON_REMAP_63        0xC8 // row address 63 is mapped to COM0


//I2Cに配列を転送する
void write_s(uint8_t *str1, uint8_t len1) {

  Wire.beginTransmission(  0x3c  );

  for (int ii = 0; ii < len1; ii++) {

    //一文字出力
    Wire.write(*str1 ++);

  }//for

  Wire.endTransmission();

}//write_s


//セットページアドレス
void setPageAddress(uint8_t start, uint8_t end)
{
  uint8_t databytes[6] = {COMMAND_MODE, SET_PAGE_ADDRESS, COMMAND_MODE, start, COMMAND_MODE, end};
  write_s(databytes, 6);
}//setPageAddress


//セットカラムアクセス
void setColumnAddress(uint8_t start, uint8_t end)
{
  uint8_t databytes[6] = {COMMAND_MODE, SET_COLUMN_ADDRESS, COMMAND_MODE, start, COMMAND_MODE, end};
  write_s(databytes, 6);
}//setColumnAddress


//セットメモリーアドレシングモード
void setMemoryAddressingMode()
{
  uint8_t databytes[4] = {COMMAND_MODE, SET_MEMORY_ADDRESSING_MODE, COMMAND_MODE, HORIZONTAL_ADDRESSING_MODE};
  write_s(databytes, 4);
}//setMemoryAddressingMode


//セットディスプレー Flip
void setDisplayFlip(int left, int down)
{
  if ( left == 0) {
    // column address   0 is mapped to SEG0 (Reset)
    //_sendCommand(SET_SEGMENT_REMAP_0);
    uint8_t databytes[2] = {COMMAND_MODE, SET_SEGMENT_REMAP_0};
    write_s(databytes, 2);
  } else {
    // column address 127 is mapped to SEG0
    //_sendCommand(SET_SEGMENT_REMAP_127);
    uint8_t databytes[2] = {COMMAND_MODE, SET_SEGMENT_REMAP_127};
    write_s(databytes, 2);
  }//end if

  if ( down == 0) {
    // Reset mode
    //_sendCommand(SET_COMMON_REMAP_0);
    uint8_t databytes[2] = {COMMAND_MODE, SET_COMMON_REMAP_0};
    write_s(databytes, 2);
  } else {
    // Flip Up/Down (Need to rewrite display before H effect shows)
    //_sendCommand(SET_COMMON_REMAP_63);
    uint8_t databytes[2] = {COMMAND_MODE, SET_COMMON_REMAP_63};
    write_s(databytes, 2);
  }//end if
}//setDisplayFlip


//8バイト分まとめて出力
void put8byte(uint8_t qq) {

  //データの配列の定義
  static uint8_t databytes[9] = {DATA_MODE, 0, 0, 0, 0, 0, 0, 0, 0};

  static char byte_count = 1;

  databytes[byte_count++] = qq;

  if ( byte_count > 8 ) { //8バイト分、貯まったら
    write_s(databytes, 9);
    byte_count = 1;
  } // endif

} //put8byte


//ドットを打つ
void Dot(int x, int y, uint8_t c) {

  static uint8_t qq = 0; //一時

  qq = qq | ( c << (y & 0x07));

  if ( (y & 0x07) == 7 ) { //1バイト分、貯まったら
    put8byte(qq); qq = 0; //SSD1306に出力
  }//end if

}//DOt


//キャラクターRAMの内容
char y_[] = {
  //1    2    3   4   5   6   7   8

    1 ,  2 ,  3 , 4 , 5 , 6 , 7 , 8  , //1
    0 ,  0 ,  0 , 0 , 0 , 0 , 0 , 0  , //2 
    0 ,  0 ,  0 , 0 , 0 , 0 , 0 , 0  , //3
    0 ,  0 ,  0 , 0 , 0 , 0 , 0 , 0  , //4

    0 ,  0 ,  0 , 0 , 0 , 0 , 0 , 0  , //5
    0 ,  0 ,  0 , 0 , 0 , 0 , 0 , 0  , //6 
    0 ,  0 ,  0 , 0 , 0 , 0 , 0 , 0  , //7
    0 ,  0 ,  0 , 0 , 0 , 0 , 0 , 0  , //8

    0 ,  0 ,  0 , 0 , 0 , 0 , 0 , 0  , //9
    0 ,  0 ,  0 , 0 , 0 , 0 , 0 , 0  , //A 
    0 ,  0 ,  0 , 0 , 0 , 0 , 0 , 0  , //B
    0 ,  0 ,  0 , 0 , 0 , 0 , 0 , 0  , //C

    0 ,  0 ,  0 , 0 , 0 , 0 , 0 , 0  , //D
    0 ,  0 ,  0 , 0 , 0 , 0 , 0 , 0  , //E 
    0 ,  0 ,  0 , 0 , 0 , 0 , 0 , 0  , //F
   56 , 57 , 58 ,59 ,60 ,61 ,62 ,63    //G

};


//再表示
void display(void) {

  char y; char x; char y1; char y2;

  //範囲の設定 (OLED内部のx,yカウンターを初期化してホームポジション0,0に)
  setPageAddress(0, MAX_PAGE);  // all pages
  setColumnAddress(0, MAX_COL); // all columns

  for (int ii = 0; ii < 8192; ii++) {

    //SSD1306のバッファーの配置順のxとyを求める
    //y =  ((ii & 0b0001110000000000 ) >> 7) + (ii & 0b0111);
    y2   =  (char)(ii       &     0b0111);
    if ( y2 == 0 ) {
      y1 =  (char)(ii >> 7) & 0b00111000;
      //x =   (ii & 0b0000001111111000)  >> 3;
      x = (char)(ii >> 3) & 0b1111111;
    }//endif y2
    y = y1 | y2;

    //↓開始　グラフのメイン処理 x=0...127 , y=0...64 
    if (y < (y_[x]+1)) {
      Dot(x, y, 1); //白のドットを打つ
    } else {
      Dot(x, y, 0); //黒のドットを打つ
    }//end if
    //↑終了

  }//for ii

}//display


//SSD1306の初期化
void display_begin(void) {

  //I2Cの初期化
  Wire.begin(); //M5Stamp S3
  delay(200);
  //Wire.setClock(2000000); //速度の変更 (I2C高速化 2Mhz) 493ms
  //Wire.setClock(1000000); //速度の変更 (I2C高速化 1Mhz) 125ms
  Wire.setClock(800000);  //速度の変更 (I2C高速化 800khz) 125ms
  //Wire.setClock(400000);  //速度の変更 (I2C高速化 400khz) 137ms
  //Wire.setClock(200000);  //速度の変更 (I2C高速化 200khz) 166ms
  //Wire.setClock(100000);  //速度の変更 (I2C高速化 100khz) 226ms
  delay(200);

  //SSD1306の初期化スペル(魔法)
  //0x80,0x8D,0x80,0x14,0x80,0xAF
  write_s( (uint8_t*) "\200\215\200\024\200\257", 6);
  delay(100);

  //セットメモリーアドレシングモード (画面の終端に来たら画面の先頭に)
  setMemoryAddressingMode();

  //セットディスプレー Flip (画面の向きを変える)
  setDisplayFlip(1, 0);

}//display_begin


//初期化
void setup() {

  // initialize serial communication at 9600 bits per second:
  Serial.begin(9600);
  Serial.println();
  Serial.println("START");
  Serial.println();

  // initialize digital pin LED_BUILTIN as an output.
  pinMode(LED_BUILTIN, OUTPUT);

  //SSD1306の初期化
  display_begin();

}//setup


//メインループ
// the loop function runs over and over again forever
void loop() {

/*
  //規則性のあるデータ
  unsigned char gg[] = {
    //0   1   2   3   4   5   6   7
    127, 255, 255, 127, 127,  0,  0, 127, //1
    127, 255, 255, 127, 127,  0,  0, 127, //2
    127, 255, 255, 127, 127,  0,  0, 127, //3
    127, 255, 255, 127, 127,  0,  0, 127, //4

    127, 255, 255, 127, 127,  0,  0, 127, //5
    127, 255, 255, 127, 127,  0,  0, 127, //6
    127, 255, 255, 127, 127,  0,  0, 127, //7
    127, 255, 255, 127, 127,  0,  0, 127, //8
  };
*/

  //ランダムウオーク debug
  static int L1 = 0;

  char i; //ループカウンター

  //入力
  summ = 0;
  for (i = 0; i < 64; i++) {

    //ランダムウオーク debug
    L1 = L1 + ( random(16) - 8 );
    if ( L1 < 0 ) {L1 = 0;} else if ( L1 > 255 ) {L1 = 255;}
    //L1 = 128;
    
    //microseconds = micros();

    //data[i] = ((analogRead(A3)) >> 2) - 128;                        // Fitting analogRead data (range:0 - 1023) to int8_t array (range:-128 - 127)
    //data_a[i] = gg[i] - 128; //規則性のあるデータ debug
    data_a[i] = L1 - 128;    //ランダムウオーク debug
    summ += data_a[i];
    //while (micros() < (microseconds + sampling_period_us)) {        // Timing out uC ADC to fulfill sampling frequency requirement
    //}
  }//for

#define DIV64(aa) (aa<0?(aa+63)>>6:aa>>6)

  //データ列から平均を引く
  // Eliminating remaining DC component (produces usable data in FFT bin #0, which is usually swamped by DC bias)
  //avg = summ / 64;
  avg = DIV64(summ);
  
  for (i = 0; i < 64; i++) {
    data_a[i] -= avg;
  }//for

  //フーリエ変換
  fix_fftr(data_a, 6, 0);                             // Performing real FFT

  char data_b[32], a, t; //一時ワーク
  for (i = 0; i < 32; i++) {

    //偶数の時は、前半 奇数の時は、後半
    a = data_a[((i << 4) & 0x10) | (i >> 1)];

    //0より小さい時は、0でそれ以外は、8倍し1/2する。(4倍)
    //63以上は、63
    if ( a < 0 ) {

      data_b[i] = 0;

    } else if (a > 15) { //a=16を4倍すると64になる為

      data_b[i] = 63;

    } else {

      data_b[i] = a * 4;

    }//if

  }//for


//SSD1306への出力　開始　↓

  for(i=0;i < 32;i++){
    
    a = data_b[i];

    t = i << 2; //x4
    y_[t++]=a;
    y_[t++]=a;
    y_[t++]=-1;
    y_[t++]=-1;

  }//for ii

  display(); //再表示

  //delay(1000); //1秒待つ debug

//SSD1306への出力　終了　↑


  digitalWrite(LED_BUILTIN, HIGH);   // turn the LED on (HIGH is the voltage level)
  delay(100);                       // wait for a second
  digitalWrite(LED_BUILTIN, LOW);    // turn the LED off by making the voltage LOW
  delay(100);  // wait for a second

  //if(data_a[2] > 16) {while(1){}} //debug
  
}//loop



```

fix_fft.cpp


```c++

#include "fix_fft.h"
/* #include <WProgram.h> */

/* fix_fft.c - Fixed-point in-place Fast Fourier Transform  */
/*
  All data are fixed-point short integers, in which -32768
  to +32768 represent -1.0 to +1.0 respectively. Integer
  arithmetic is used for speed, instead of the more natural
  floating-point.

  For the forward FFT (time -> freq), fixed scaling is
  performed to prevent arithmetic overflow, and to map a 0dB
  sine/cosine wave (i.e. amplitude = 32767) to two -6dB freq
  coefficients. The return value is always 0.

  For the inverse FFT (freq -> time), fixed scaling cannot be
  done, as two 0dB coefficients would sum to a peak amplitude
  of 64K, overflowing the 32k range of the fixed-point integers.
  Thus, the fix_fft() routine performs variable scaling, and
  returns a value which is the number of bits LEFT by which
  the output must be shifted to get the actual amplitude
  (i.e. if fix_fft() returns 3, each value of fr[] and fi[]
  must be multiplied by 8 (2**3) for proper scaling.
  Clearly, this cannot be done within fixed-point short
  integers. In practice, if the result is to be used as a
  filter, the scale_shift can usually be ignored, as the
  result will be approximately correctly normalized as is.

  Written by:  Tom Roberts  11/8/89
  Made portable:  Malcolm Slaney 12/15/94 malcolm@interval.com
  Enhanced:  Dimitrios P. Bouras  14 Jun 2006 dbouras@ieee.org
  Modified for 8bit values David Keller  10.10.2010
*/

/*
  FIX_MPY() - fixed-point multiplication & scaling.
  Substitute inline assembly for hardware-specific
  optimization suited to a particluar DSP processor.
  Scaling ensures that result remains 16-bit.
*/
inline int8_t FIX_MPY(int8_t a, int8_t b)
{

  //Serial.println(a);
 //Serial.println(b);


    /* shift right one less bit (i.e. 15-1) */
    int16_t c = ((int16_t)a * (int16_t)b) >> 6;
    /* last bit shifted out = rounding-bit */
    b = c & 0x01;
    /* last shift + rounding bit */
    a = (c >> 1) + b;

      /*
      Serial.println(Sinewave[3]);
      Serial.println(c);
      Serial.println(a);
      while(1);*/

    return a;
}

/*
  fix_fft() - perform forward/inverse fast Fourier transform.
  fr[n],fi[n] are real and imaginary arrays, both INPUT AND
  RESULT (in-place FFT), with 0 <= n < 2**m; set inverse to
  0 for forward transform (FFT), or 1 for iFFT.
*/
int16_t fix_fft(int8_t fr[], int8_t fi[], int16_t m, int16_t inverse)
{
    int16_t mr, nn, i, j, l, k, istep, n, scale, shift;
    int8_t qr, qi, tr, ti, wr, wi;

    n = 1 << m;

    /* max FFT size = N_WAVE */
    if (n > N_WAVE)
      return -1;

    mr = 0;
    nn = n - 1;
    scale = 0;

    /* decimation in time - re-order data */
    for (m=1; m<=nn; ++m) {
      l = n;
      do {
        l >>= 1;
      } while (mr+l > nn);
      mr = (mr & (l-1)) + l;

      if (mr <= m)
        continue;
      tr = fr[m];
      fr[m] = fr[mr];
      fr[mr] = tr;
      ti = fi[m];
      fi[m] = fi[mr];
      fi[mr] = ti;
    }

    l = 1;
    k = LOG2_N_WAVE-1;
    while (l < n) {
      if (inverse) {
        /* variable scaling, depending upon data */
        shift = 0;
        for (i=0; i<n; ++i) {
            j = fr[i];
            if (j < 0)
              j = -j;
            m = fi[i];
            if (m < 0)
              m = -m;
            if (j > 16383 || m > 16383) {
              shift = 1;
              break;
            }
        }
        if (shift)
            ++scale;
      } else {
        /*
          fixed scaling, for proper normalization --
          there will be log2(n) passes, so this results
          in an overall factor of 1/n, distributed to
          maximize arithmetic accuracy.
        */
        shift = 1;
      }
      /*
        it may not be obvious, but the shift will be
        performed on each data point exactly once,
        during this pass.
      */
      istep = l << 1;
      for (m=0; m<l; ++m) {
        j = m << k;
        /* 0 <= j < N_WAVE/2 */
        wr =  pgm_read_byte_near(Sinewave + j+N_WAVE/4);

/*Serial.println("asdfasdf");
Serial.println(wr);
Serial.println(j+N_WAVE/4);
Serial.println(Sinewave[256]);

Serial.println("");*/


        wi = -pgm_read_byte_near(Sinewave + j);
        if (inverse)
            wi = -wi;
        if (shift) {
            wr >>= 1;
            wi >>= 1;
        }
        for (i=m; i<n; i+=istep) {
            j = i + l;
            tr = FIX_MPY(wr,fr[j]) - FIX_MPY(wi,fi[j]);
            ti = FIX_MPY(wr,fi[j]) + FIX_MPY(wi,fr[j]);
            qr = fr[i];
            qi = fi[i];
            if (shift) {
              qr >>= 1;
              qi >>= 1;
            }
            fr[j] = qr - tr;
            fi[j] = qi - ti;
            fr[i] = qr + tr;
            fi[i] = qi + ti;
        }
      }
      --k;
      l = istep;
    }
    return scale;
}

/*
  fix_fftr() - forward/inverse FFT on array of real numbers.
  Real FFT/iFFT using half-size complex FFT by distributing
  even/odd samples into real/imaginary arrays respectively.
  In order to save data space (i.e. to avoid two arrays, one
  for real, one for imaginary samples), we proceed in the
  following two steps: a) samples are rearranged in the real
  array so that all even samples are in places 0-(N/2-1) and
  all imaginary samples in places (N/2)-(N-1), and b) fix_fft
  is called with fr and fi pointing to index 0 and index N/2
  respectively in the original array. The above guarantees
  that fix_fft "sees" consecutive real samples as alternating
  real and imaginary samples in the complex array.
*/
int16_t fix_fftr(int8_t f[], int16_t m, int16_t inverse)
{
    int16_t i, N = 1<<(m-1), scale = 0;
    int8_t tt, *fr=f, *fi=&f[N];

    if (inverse)
      scale = fix_fft(fi, fr, m-1, inverse);
    for (i=1; i<N; i+=2) {
      tt = f[N+i-1];
      f[N+i-1] = f[i];
      f[i] = tt;
    }
    if (! inverse)
      scale = fix_fft(fi, fr, m-1, inverse);
    return scale;
}


```

fix_fft.h


```c++

#ifndef FIXFFT_H
#define FIXFFT_H

#if (defined(__AVR__))
  #include <avr/pgmspace.h>
#else
  #include <pgmspace.h>
#endif

#if ARDUINO >= 100
  #include "Arduino.h"
#else
  #include "WProgram.h" /* This is where the standard Arduino code lies */
#endif

#define N_WAVE  256    /* full length of Sinewave[] */
#define LOG2_N_WAVE 8   /* log2(N_WAVE) */

/*
  Since we only use 3/4 of N_WAVE, we define only
  this many samples, in order to conserve data space.
*/

const int8_t Sinewave[N_WAVE-N_WAVE/4] PROGMEM = {
0, 3, 6, 9, 12, 15, 18, 21,
24, 28, 31, 34, 37, 40, 43, 46,
48, 51, 54, 57, 60, 63, 65, 68,
71, 73, 76, 78, 81, 83, 85, 88,
90, 92, 94, 96, 98, 100, 102, 104,
106, 108, 109, 111, 112, 114, 115, 117,
118, 119, 120, 121, 122, 123, 124, 124,
125, 126, 126, 127, 127, 127, 127, 127,

127, 127, 127, 127, 127, 127, 126, 126,
125, 124, 124, 123, 122, 121, 120, 119,
118, 117, 115, 114, 112, 111, 109, 108,
106, 104, 102, 100, 98, 96, 94, 92,
90, 88, 85, 83, 81, 78, 76, 73,
71, 68, 65, 63, 60, 57, 54, 51,
48, 46, 43, 40, 37, 34, 31, 28,
24, 21, 18, 15, 12, 9, 6, 3,

0, -3, -6, -9, -12, -15, -18, -21,
-24, -28, -31, -34, -37, -40, -43, -46,
-48, -51, -54, -57, -60, -63, -65, -68,
-71, -73, -76, -78, -81, -83, -85, -88,
-90, -92, -94, -96, -98, -100, -102, -104,
-106, -108, -109, -111, -112, -114, -115, -117,
-118, -119, -120, -121, -122, -123, -124, -124,
-125, -126, -126, -127, -127, -127, -127, -127,

/*-127, -127, -127, -127, -127, -127, -126, -126,
-125, -124, -124, -123, -122, -121, -120, -119,
-118, -117, -115, -114, -112, -111, -109, -108,
-106, -104, -102, -100, -98, -96, -94, -92,
-90, -88, -85, -83, -81, -78, -76, -73,
-71, -68, -65, -63, -60, -57, -54, -51,
-48, -46, -43, -40, -37, -34, -31, -28,
-24, -21, -18, -15, -12, -9, -6, -3, */
};

/*
  fix_fft() - perform forward/inverse fast Fourier transform.
  fr[n],fi[n] are real and imaginary arrays, both INPUT AND
  RESULT (in-place FFT), with 0 <= n < 2**m; set inverse to
  0 for forward transform (FFT), or 1 for iFFT.
*/
int16_t fix_fft(int8_t fr[], int8_t fi[], int16_t m, int16_t inverse);

/*
  fix_fftr() - forward/inverse FFT on array of real numbers.
  Real FFT/iFFT using half-size complex FFT by distributing
  even/odd samples into real/imaginary arrays respectively.
  In order to save data space (i.e. to avoid two arrays, one
  for real, one for imaginary samples), we proceed in the
  following two steps: a) samples are rearranged in the real
  array so that all even samples are in places 0-(N/2-1) and
  all imaginary samples in places (N/2)-(N-1), and b) fix_fft
  is called with fr and fi pointing to index 0 and index N/2
  respectively in the original array. The above guarantees
  that fix_fft "sees" consecutive real samples as alternating
  real and imaginary samples in the complex array.
*/
int16_t fix_fftr(int8_t f[], int16_t m, int16_t inverse);

#endif


```