参考
x 軽量高速化、みやすさ重視
x 桁あふれを修正(オーバーフロー)
目的
フーリエ変換の速度を測る
結果
//fft_test3_UNO_1
//インクルド
#include "Arduino.h"
#include "fix_fft.h" // 8-bit FFT library modified for Arduino
//定義
//int8_t data[64], buff[32]; // used to store FFT input/output and past data
char data_a[64], buff_a[32];
//unsigned long microseconds; // used for timekeeping
int summ, avg; // used for DC bias elimination
//初期化
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);
}//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
};
Serial.println("START>>");
unsigned long time;
time = millis();
//入力
summ = 0;
for (int i = 0; i < 64; i++) {
//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;
summ += data_a[i];
//while (micros() < (microseconds + sampling_period_us)) { // Timing out uC ADC to fulfill sampling frequency requirement
//}
}//for
//データ列から平均を引く
// Eliminating remaining DC component (produces usable data in FFT bin #0, which is usually swamped by DC bias)
//avg = summ / 64;
if (summ < 0) {avg = -((-summ) >> 6);} else {avg = summ >> 6;}
for (int i = 0; i < 64; i++) {
data_a[i] -= avg;
}//for
//フーリエ変換
fix_fftr(data_a, 6, 0); // Performing real FFT
/*
//結果の表示
Serial.println(">>>");
int uu = 1;
for (int i = 0; i < 64; i++) {
//printf("%d,",data_a[i]);
Serial.print((int)data_a[i]);
Serial.print(',');
uu++;
if (uu > 8) {
uu = 1;
Serial.println();
}
}//for
Serial.println("<<<");
digitalWrite(LED_BUILTIN, HIGH); // turn the LED on (HIGH is the voltage level)
delay(1000); // wait for a second
digitalWrite(LED_BUILTIN, LOW); // turn the LED off by making the voltage LOW
delay(1000); // wait for a second
*/
/*
char data_a[64] = {
//1 2 3 4 5 6 7 8
-1, 0, 0, 0, -44, -1, -1, 0, //1
0, 0, -1, -1, 0, -1, -1, 0, //2
-1, 0, 0, 0, 44, -1, -1, 0, //3
-2, 0, -1, -1, -2, -1, -1, -2, //4
-1, -2, -2, -1, 106, 0, 0, 0, //5
-2, 0, 0, -1, -1, 0, 0, 0, //6
-1, 0, 0, -1, 18, 0, 0, 0, //7
0, 0, 0, -1, 1, 0, 0, 0 //8
};
*/
/*
char data_a[64] = {
//1 2 3 4 5 6 7 8
1, 2, 3, 4, 5, 6, 7, 8, //1
9, 10, 11, 12, 13, 14, 15, 16, //2
0, 0, 0, 0, 0, 0, 0, 0, //3
0, 0, 0, 0, 0, 0, 0, 0, //4
-1, -2, -2, -1, 106, 0, 0, 0, //5
-2, 0, 0, -1, -1, 0, 0, 0, //6
-1, 0, 0, -1, 18, 0, 0, 0, //7
0, 0, 0, -1, 1, 0, 0, 0 //8
};
*/
/*
char data_a[64] = {
//1 2 3 4 5 6 7 8
0, 0, 0, 0, 0, 0, 0, 0, //1
0, 0, 0, 0, 0, 0, 0, 0, //2
1, 2, 3, 4, 5, 6, 7, 8, //3
9, 10, 11, 12, 13, 14, 15, 16, //4
-1, -2, -2, -1, 106, 0, 0, 0, //5
-2, 0, 0, -1, -1, 0, 0, 0, //6
-1, 0, 0, -1, 18, 0, 0, 0, //7
0, 0, 0, -1, 1, 0, 0, 0 //8
};
*/
//int8_t reordered[32];
char reordered[32], a;
for (int i = 0; i < 32; i++) {
//偶数の時は、前半 奇数の時は、後半
a = data_a[((i << 4) & 0x10) | (i >> 1)];
//0より小さい時は、0でそれ以外は、8倍し1/2する。(4倍)
//63以上は、63
if ( a < 0 ) {
reordered[i] = 0;
} else if (a > 15) { //a=16を4倍すると64になる為
reordered[i] = 63;
} else {
reordered[i] = a * 4;
}//if
}//for
//データにコピーする
for (int i = 0; i < 32; i++) {
data_a[i] = reordered[i];
}//for
time = millis() - time;
Serial.print("millis=");
Serial.println((int)time);
//結果の表示
//Serial.println(">>>");
int uu = 1;
for (int i = 0; i < 32; i++) {
//printf("%d,",data_a[i]);
Serial.print((int)data_a[i]);
Serial.print(',');
uu++;
if (uu > 8) {
uu = 1;
//printf("\n");
Serial.println();
}//if
}//for
//Serial.println("<<<");
Serial.println("<<END");
digitalWrite(LED_BUILTIN, HIGH); // turn the LED on (HIGH is the voltage level)
delay(1000); // wait for a second
digitalWrite(LED_BUILTIN, LOW); // turn the LED off by making the voltage LOW
delay(1000); // wait for a second
}//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