1
0

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

Last updated at Posted at 2024-08-28

x 過去ログをみよ

結果

o_coq349.jpg

プログラム




//
//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


#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


#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


1
0
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
1
0