/**
*
* Copyright (c) 2006 Henrik Sundt
*
* This code is free software; you can redistribute it and/or
* modify it under the terms of the GNU Lesser General Public
* License as published by the Free Software Foundation; either
* version 2.1 of the License, or (at your option) any later version.
*
* This code is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
* Lesser General Public License for more details.
*
* More information about GNU Lesser General Public License
* can be obtained from the
* Free Software Foundation, Inc., 59 Temple Place, Suite 330,
* Boston, MA 02111-1307 USA
*
* @author Henrik Sundt
*/
// Realtime frequency analysis and display from microphone input for mobile phone.
// Analysis while recording of audio, with 1 second delay.
// User can select Goertzel or FFT algorithm alternatively.
// Both give same result, but calculates in different ways.
// This code is written for Mobile Processing ( http://mobile.processing.org )
// and Java j2me
// Written by Henrik Sundt (hsundt at notam02.no)
// October 2006
// for Trolley Singers Project (see http://silver.avu.cz/trolley_singers )
String version="v0.4r";
import javax.microedition.media.*;
import javax.microedition.media.control.*;
import java.io.*;
int amode=1; //amode==1 gives Goertzel analysis, amode==2 gives FFT analysis
long w,h;
PFont fontg,fontw,fontstd;
boolean lredraw,mredraw;
int[] Wre = new int[8];
int[] Wim_n = new int[8];
int[] result = new int[8];
String infoline,errorline,memline;
int count=0;
public class myByteArrayOutputStream extends ByteArrayOutputStream {
public myByteArrayOutputStream() {
super();
}
public myByteArrayOutputStream(int size) {
super(size);
}
public int mycount() {
return count;
}
public byte[] getBuffer() {
return buf;
}
public void bufferCopyTo(byte[] inbuf,int start,int len) {
int a;
if(start+len>=count) len=(count-start)-1;
for(a=0; a < len; a++)
inbuf[a] = buf[start+a];
}
// "int"-version of the above function:
public void bufferCopyTo(int[] inbuf,int start,int len) {
int a;
if(start+len>=count) len=(count-start)-1;
for(a=0; a < len; a++)
inbuf[a] = (int)(buf[start+a]);
}
}
public myByteArrayOutputStream[] output = new myByteArrayOutputStream[2]; // record buffer and play buffer, to be switched
Player player;
RecordControl rc;
int recorder_starttime=0;
int buf=0;
void setup() {
errorline= "everything running well";
infoline= "everything running well";
memline="";
framerate(99);
lredraw = false;
mredraw=false;
fontg = loadFont("verdana.mvlw",color(100,255,100));
fontw = loadFont("verdana.mvlw",color(200,200,200));
fontstd = loadFont(FACE_MONOSPACE, STYLE_PLAIN, SIZE_SMALL);
init_bars();
init_fft();
// Program assumes 8000 Hz samplerate
// Formulas used for computing frequency: Change k for other frequencies.
// 12 bit fixed point calculations.
//
// Wre=cos(2*PI*k/N)*(1 << 12) N=8000 (samplerate) k= # Hz
// Wim_n=sin(2*PI*k/N)*(1<<12)
// (1<<12) = 4096
//
Wre[0] = 3983; // 300 Hz
Wim_n[0] = 956;
Wre[1] = 3896; // 400 Hz
Wim_n[1] = 1266;
Wre[2] = 3784; // 500 Hz
Wim_n[2] = 1567;
Wre[3] = 3650; // 600 Hz
Wim_n[3] = 1860;
Wre[4] = 3492; // 700 Hz
Wim_n[4] = 2140;
Wre[5] = 3314; // 800 Hz
Wim_n[5] = 2408;
Wre[6] = 3115; // 900 Hz
Wim_n[6] = 2660;
Wre[7] = 2896; // 1000 Hz
Wim_n[7] = 2896;
for(int a=0; a<8; a++) { // adjust fix point to 10 instead of 12 bits (used later)
Wre[a] /= (1<<2);
Wim_n[a] /= (1<<2);
}
w = width;
h = height;
softkey("");
try {
output[0] = new myByteArrayOutputStream(4*8000); // buffer
output[1] = new myByteArrayOutputStream(4*8000); // buffer
} catch(Exception e) {
e.printStackTrace();
errorline="init failed";
}
mredraw=true;
}
//------------------------------------------------
int ctime=0;
Thread mythread;
boolean newthread=true;
void draw()
{
// this code prevents second thread generated by keypress to run (bug in some versions of Processing)
Thread thisthread=Thread.currentThread();
if(newthread) {
mythread=thisthread;
newthread=false;
}
else if(thisthread != mythread) {
try {
Thread.currentThread().sleep(10000); // unneeded thread sleeps for 10 secs
} catch(InterruptedException ie)
{}
return;
}
// -- end ----------------------------------------------------
// Start of the actual draw() :
int pos;
int mtime=0;
int currtime=millis();
if(lredraw) {
if(ctime==0) ctime=currtime;
mtime=currtime-recorder_starttime; // current relative recorder position, in milliseconds
pos=(int)((mtime*8)); // convert to sample position at 8000 rate
if(amode==1)
goertzel((buf*-1)+1,pos); // do analysis on the opposite buffer, same position as record position
else if(amode==2)
fft_analysis((buf*-1)+1,pos);
printresult(pos);
}
else if(mredraw) {
menu();
mredraw=false;
}
if(lredraw && (currtime-ctime > 1000)) { // restart recording and switch record buffer
try {
close_record();
buf=(buf*-1)+1; // switch buffer from 0 to 1, or 1 to 0
start_record(buf);
memline="free mem: "+str(freeMemory());
} catch(Exception e) {
e.printStackTrace();
errorline="recorder restart failed";
}
ctime=millis();
}
try {
Thread.currentThread().sleep(50); // delay between analysises // must be >=40 on N70 phone
} catch(InterruptedException ie) {
errorline= "sleep in draw() failed";
}
}
// destructor
void destroy()
{
try {
rc.commit();
} catch(Exception e) {
e.printStackTrace();
}
try {
player.close(); // needed because otherwise program will hang on phone after exit
} catch(Exception e) {
e.printStackTrace();
}
}
//------------------------------------------------
void menu(){
background(0);
stroke(0,80,0);
fill(0,40,0);
rect(-5,14,width+10,30);
textAlign(LEFT);
textFont(fontg);
text("Goertzel/fft realtime "+version,20,30);
text("Press 1 to start goertzel",20,60);
text("Press 2 to start fft",20,72);
text("with realtime analysis",20,84);
}
//------------------------------------------------
void playbacksound(int milliseconds, byte[] wav) throws IOException, MediaException, InterruptedException
{
Player audio;
audio = Manager.createPlayer(new ByteArrayInputStream(wav),"audio/x-wav");
audio.prefetch();
audio.start();
Thread.currentThread().sleep(milliseconds);
audio.close();
}
void start_record(int bb)
{
try {
player=Manager.createPlayer("capture://audio?encoding=pcm&bits=8&rate=8000");
player.realize();
rc = (RecordControl)player.getControl("RecordControl");
} catch(Exception e) {
e.printStackTrace();
errorline="rec failed (2.5)";
}
try {
output[bb].reset();
} catch(Exception e) {
e.printStackTrace();
errorline="rec failed (3)";
}
try {
rc.setRecordStream(output[bb]);
} catch(Exception e) {
e.printStackTrace();
errorline="rec failed (4)";
}
try {
rc.startRecord();
} catch(Exception e) {
e.printStackTrace();
errorline="rec failed (5)";
}
try {
player.start();
recorder_starttime=millis();
} catch(Exception e) {
e.printStackTrace();
errorline="rec failed (6)";
}
}
void close_record()
{
try {
rc.commit();
} catch(Exception e) {
e.printStackTrace();
errorline="rec failed (8)";
}
try {
player.stop();
} catch(Exception e) {
e.printStackTrace();
errorline="rec failed (9)";
}
try {
player.close();
} catch(Exception e) {
e.printStackTrace();
errorline="rec failed (10)";
}
}
// GOERTZEL ALGORITHM:
int v,v1,v2;
int g_iter=250; // number of iterations in Goertzel
byte[] signal = new byte[g_iter];
void goertzel(int bb, int pos)
{
int x,f,sample;
int yre,yim;
int sz=0;
sz=output[bb].size();
if(pos < 0 || 44+pos+g_iter>=sz) return;
infoline="bb="+bb+" r="+str(sz)+" c="+output[bb].mycount();
output[bb].bufferCopyTo(signal,44+pos,g_iter); // skipping header in .wav-data (temporary solution)
for(f=0; f < 8; f++) {
v1=0; v2=0;
for(x=0; x < g_iter; x++) {
// read a sample from .wav-data
// (efficiency in loop could be improved ...)
sample=(int)(((signal[x]))); // 8 bit samples
v=((Wre[f]*v1)>>(10-1))-v2+(sample);
v2=v1;
v1=v;
}
yre=v-((Wre[f]*v2)>>10);
yim=(Wim_n[f]*v2)>>10;
//result[f] = mylog((yre*yre+yim*yim)/g_iter); // result in log
result[f] = (((yre*yre+yim*yim)>>12)/g_iter); // result in linear
}
}
// PLOTTING BARS:
String[] bar = new String[21];
void init_bars()
{
bar[0]="";
bar[1]="*";
bar[2]="**";
bar[3]="***";
bar[4]="****";
bar[5]="*****";
bar[6]="******";
bar[7]="*******";
bar[8]="********";
bar[9]="*********";
bar[10]="**********";
bar[11]="***********";
bar[12]="************";
bar[13]="*************";
bar[14]="**************";
bar[15]="***************";
bar[16]="****************";
bar[17]="*****************";
bar[18]="******************";
bar[19]="*******************";
bar[20]="********************";
}
int[] bresult = new int[8];
void printresult(int pos)
{
String st;
for(int x=0; x < 8; x++) { // limit the bar size to 20
bresult[x]=(result[x]);
if(bresult[x]>20) bresult[x]=20;
}
background(0);
stroke(255);
fill(255);
textAlign(LEFT);
textFont(fontstd);
st=" 300 "+bar[bresult[0]]; text(st,0,24);
st=" 400 "+bar[bresult[1]]; text(st,0,24+12*1);
st=" 500 "+bar[bresult[2]]; text(st,0,24+12*2);
st=" 600 "+bar[bresult[3]]; text(st,0,24+12*3);
st=" 700 "+bar[bresult[4]]; text(st,0,24+12*4);
st=" 800 "+bar[bresult[5]]; text(st,0,24+12*5);
st=" 900 "+bar[bresult[6]]; text(st,0,24+12*6);
st="1000 "+bar[bresult[7]]; text(st,0,24+12*7);
text(infoline,20,108+12);
st="count "+str(count)+" pos "+str(pos);
text(st,20,108+12*2);
count++;
text(memline,20,108+12*3);
if(amode==1)
text("goert",140,24);
else if(amode==2)
text("fft",140,24);
}
void keyReleased()
{
if(!lredraw) { start_record(buf); lredraw=true; }
if(key=='1') amode=1;
if(key=='2') amode=2;
}
////////////////////// FFT //////////////////////////////////////////
int M=9;
int N=1<=sz) return;
infoline="bb="+bb+" r="+str(sz)+" c="+output[bb].mycount();
output[bb].bufferCopyTo(real,44+pos,N);
fix_fft(real, imag, M, false);
// 300 hz = bin 19.2
// 400 hz = bin 25.6
// 500 hz = bin 32
// 600 hz = bin 38.4
// 700 hz = bin 44.8
// 800 hz = bin 51.2
// 900 hz = bin 57.6
// 1000hz = bin 64
//
// approximations:
/* logarithmic display:
result[0] = mylog(real[19]*real[19]+imag[19]*imag[19]);
result[1] = mylog(real[26]*real[26]+imag[26]*imag[26]);
result[2] = mylog(real[32]*real[32]+imag[32]*imag[32]);
result[3] = mylog(real[38]*real[38]+imag[38]*imag[38]);
result[4] = mylog(real[45]*real[45]+imag[45]*imag[45]);
result[5] = mylog(real[51]*real[51]+imag[51]*imag[51]);
result[6] = mylog(real[57]*real[57]+imag[57]*imag[57]);
result[7] = mylog(real[64]*real[64]+imag[64]*imag[64]);
*/
// linear display of power spectrum
result[0] = (real[19]*real[19]+imag[19]*imag[19])>>4;
result[1] = (real[26]*real[26]+imag[26]*imag[26])>>4;
result[2] = (real[32]*real[32]+imag[32]*imag[32])>>4;
result[3] = (real[38]*real[38]+imag[38]*imag[38])>>4;
result[4] = (real[45]*real[45]+imag[45]*imag[45])>>4;
result[5] = (real[51]*real[51]+imag[51]*imag[51])>>4;
result[6] = (real[57]*real[57]+imag[57]*imag[57])>>4;
result[7] = (real[64]*real[64]+imag[64]*imag[64])>>4;
}
int mylog(int a) // an approximate and quick 2-logarithm calculation
{
int logish=0;
while(a>0) {
a/=2;
logish++;
}
return logish;
}
/* fix_fft.c - Fixed-point Fast Fourier Transform */
/*
fix_fft() perform FFT or inverse FFT
window() applies a Hanning window to the (time) input
fix_loud() calculates the loudness of the signal, for
each freq point. Result is an integer array,
units are dB (values will be negative).
iscale() scale an integer value by (numer/denom).
fix_mpy() perform fixed-point multiplication.
Sinewave[1024] sinewave normalized to 32767 (= 1.0).
Loudampl[100] Amplitudes for lopudnesses from 0 to -99 dB.
Low_pass Low-pass filter, cutoff at sample_freq / 4.
All data are fixed-point short integers, in which
-32768 to +32768 represent -1.0 to +1.0. 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 one in the lower half is reported as 0dB
by fix_loud(). 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 the 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.
TURBO C, any memory model; uses inline assembly for speed
and for carefully-scaled arithmetic.
Written by: Tom Roberts 11/8/89
Made portable: Malcolm Slaney 12/15/94 malcolm at interval.com
Timing on a Macintosh PowerBook 180.... (using Symantec C6.0)
fix_fft (1024 points) 8 ticks
fft (1024 points - Using SANE) 112 Ticks
fft (1024 points - Using FPU) 11
*/
int[] Sinewave = {
0, 201, 402, 603, 804, 1005, 1206, 1406,
1607, 1808, 2009, 2209, 2410, 2610, 2811, 3011,
3211, 3411, 3611, 3811, 4011, 4210, 4409, 4608,
4807, 5006, 5205, 5403, 5601, 5799, 5997, 6195,
6392, 6589, 6786, 6982, 7179, 7375, 7571, 7766,
7961, 8156, 8351, 8545, 8739, 8932, 9126, 9319,
9511, 9703, 9895, 10087, 10278, 10469, 10659, 10849,
11038, 11227, 11416, 11604, 11792, 11980, 12166, 12353,
12539, 12724, 12909, 13094, 13278, 13462, 13645, 13827,
14009, 14191, 14372, 14552, 14732, 14911, 15090, 15268,
15446, 15623, 15799, 15975, 16150, 16325, 16499, 16672,
16845, 17017, 17189, 17360, 17530, 17699, 17868, 18036,
18204, 18371, 18537, 18702, 18867, 19031, 19194, 19357,
19519, 19680, 19840, 20000, 20159, 20317, 20474, 20631,
20787, 20942, 21096, 21249, 21402, 21554, 21705, 21855,
22004, 22153, 22301, 22448, 22594, 22739, 22883, 23027,
23169, 23311, 23452, 23592, 23731, 23869, 24006, 24143,
24278, 24413, 24546, 24679, 24811, 24942, 25072, 25201,
25329, 25456, 25582, 25707, 25831, 25954, 26077, 26198,
26318, 26437, 26556, 26673, 26789, 26905, 27019, 27132,
27244, 27355, 27466, 27575, 27683, 27790, 27896, 28001,
28105, 28208, 28309, 28410, 28510, 28608, 28706, 28802,
28897, 28992, 29085, 29177, 29268, 29358, 29446, 29534,
29621, 29706, 29790, 29873, 29955, 30036, 30116, 30195,
30272, 30349, 30424, 30498, 30571, 30643, 30713, 30783,
30851, 30918, 30984, 31049,
31113, 31175, 31236, 31297,
31356, 31413, 31470, 31525, 31580, 31633, 31684, 31735,
31785, 31833, 31880, 31926, 31970, 32014, 32056, 32097,
32137, 32176, 32213, 32249, 32284, 32318, 32350, 32382,
32412, 32441, 32468, 32495, 32520, 32544, 32567, 32588,
32609, 32628, 32646, 32662, 32678, 32692, 32705, 32717,
32727, 32736, 32744, 32751, 32757, 32761, 32764, 32766,
32767, 32766, 32764, 32761, 32757, 32751, 32744, 32736,
32727, 32717, 32705, 32692, 32678, 32662, 32646, 32628,
32609, 32588, 32567, 32544, 32520, 32495, 32468, 32441,
32412, 32382, 32350, 32318, 32284, 32249, 32213, 32176,
32137, 32097, 32056, 32014, 31970, 31926, 31880, 31833,
31785, 31735, 31684, 31633, 31580, 31525, 31470, 31413,
31356, 31297, 31236, 31175, 31113, 31049, 30984, 30918,
30851, 30783, 30713, 30643, 30571, 30498, 30424, 30349,
30272, 30195, 30116, 30036, 29955, 29873, 29790, 29706,
29621, 29534, 29446, 29358, 29268, 29177, 29085, 28992,
28897, 28802, 28706, 28608, 28510, 28410, 28309, 28208,
28105, 28001, 27896, 27790, 27683, 27575, 27466, 27355,
27244, 27132, 27019, 26905, 26789, 26673, 26556, 26437,
26318, 26198, 26077, 25954, 25831, 25707, 25582, 25456,
25329, 25201, 25072, 24942, 24811, 24679, 24546, 24413,
24278, 24143, 24006, 23869, 23731, 23592, 23452, 23311,
23169, 23027, 22883, 22739, 22594, 22448, 22301, 22153,
22004, 21855, 21705, 21554, 21402, 21249, 21096, 20942,
20787, 20631, 20474, 20317, 20159, 20000, 19840, 19680,
19519, 19357, 19194, 19031, 18867, 18702, 18537, 18371,
18204, 18036, 17868, 17699, 17530, 17360, 17189, 17017,
16845, 16672, 16499, 16325, 16150, 15975, 15799, 15623,
15446, 15268, 15090, 14911, 14732, 14552, 14372, 14191,
14009, 13827, 13645, 13462, 13278, 13094, 12909, 12724,
12539, 12353, 12166, 11980, 11792, 11604, 11416, 11227,
11038, 10849, 10659, 10469, 10278, 10087, 9895, 9703,
9511, 9319, 9126, 8932, 8739, 8545, 8351, 8156,
7961, 7766, 7571, 7375, 7179, 6982, 6786, 6589,
6392, 6195, 5997, 5799, 5601, 5403, 5205, 5006,
4807, 4608, 4409, 4210, 4011, 3811, 3611, 3411,
3211, 3011, 2811, 2610, 2410, 2209, 2009, 1808,
1607, 1406, 1206, 1005, 804, 603, 402, 201,
0, -201, -402, -603, -804, -1005, -1206, -1406,
-1607, -1808, -2009, -2209, -2410, -2610, -2811, -3011,
-3211, -3411, -3611, -3811, -4011, -4210, -4409, -4608,
-4807, -5006, -5205, -5403, -5601, -5799, -5997, -6195,
-6392, -6589, -6786, -6982, -7179, -7375, -7571, -7766,
-7961, -8156, -8351, -8545, -8739, -8932, -9126, -9319,
-9511, -9703, -9895, -10087, -10278, -10469, -10659, -10849,
-11038, -11227, -11416, -11604, -11792, -11980, -12166, -12353,
-12539, -12724, -12909, -13094, -13278, -13462, -13645, -13827,
-14009, -14191, -14372, -14552, -14732, -14911, -15090, -15268,
-15446, -15623, -15799, -15975, -16150, -16325, -16499, -16672,
-16845, -17017, -17189, -17360, -17530, -17699, -17868, -18036,
-18204, -18371, -18537, -18702, -18867, -19031, -19194, -19357,
-19519, -19680, -19840, -20000, -20159, -20317, -20474, -20631,
-20787, -20942, -21096, -21249, -21402, -21554, -21705, -21855,
-22004, -22153, -22301, -22448, -22594, -22739, -22883, -23027,
-23169, -23311, -23452, -23592, -23731, -23869, -24006, -24143,
-24278, -24413, -24546, -24679, -24811, -24942, -25072, -25201,
-25329, -25456, -25582, -25707, -25831, -25954, -26077, -26198,
-26318, -26437, -26556, -26673, -26789, -26905, -27019, -27132,
-27244, -27355, -27466, -27575, -27683, -27790, -27896, -28001,
-28105, -28208, -28309, -28410, -28510, -28608, -28706, -28802,
-28897, -28992, -29085, -29177, -29268, -29358, -29446, -29534,
-29621, -29706, -29790, -29873, -29955, -30036, -30116, -30195,
-30272, -30349, -30424, -30498, -30571, -30643, -30713, -30783,
-30851, -30918, -30984, -31049, -31113, -31175, -31236, -31297,
-31356, -31413, -31470, -31525, -31580, -31633, -31684, -31735,
-31785, -31833, -31880, -31926, -31970, -32014, -32056, -32097,
-32137, -32176, -32213, -32249, -32284, -32318, -32350, -32382,
-32412, -32441, -32468, -32495, -32520, -32544, -32567, -32588,
-32609, -32628, -32646, -32662, -32678, -32692, -32705, -32717,
-32727, -32736, -32744, -32751, -32757, -32761, -32764, -32766,
-32767, -32766, -32764, -32761, -32757, -32751, -32744, -32736,
-32727, -32717, -32705, -32692, -32678, -32662, -32646, -32628,
-32609, -32588, -32567, -32544, -32520, -32495, -32468, -32441,
-32412, -32382, -32350, -32318, -32284, -32249, -32213, -32176,
-32137, -32097, -32056, -32014, -31970, -31926, -31880, -31833,
-31785, -31735, -31684, -31633, -31580, -31525, -31470, -31413,
-31356, -31297, -31236, -31175, -31113, -31049, -30984, -30918,
-30851, -30783, -30713, -30643, -30571, -30498, -30424, -30349,
-30272, -30195, -30116, -30036, -29955, -29873, -29790, -29706,
-29621, -29534, -29446, -29358, -29268, -29177, -29085, -28992,
-28897, -28802, -28706, -28608, -28510, -28410, -28309, -28208,
-28105, -28001, -27896, -27790, -27683, -27575, -27466, -27355,
-27244, -27132, -27019, -26905, -26789, -26673, -26556, -26437,
-26318, -26198, -26077, -25954, -25831, -25707, -25582, -25456,
-25329, -25201, -25072, -24942, -24811, -24679, -24546, -24413,
-24278, -24143, -24006, -23869, -23731, -23592, -23452, -23311,
-23169, -23027, -22883, -22739, -22594, -22448, -22301, -22153,
-22004, -21855, -21705, -21554, -21402, -21249, -21096, -20942,
-20787, -20631, -20474, -20317, -20159, -20000, -19840, -19680,
-19519, -19357, -19194, -19031, -18867, -18702, -18537, -18371,
-18204, -18036, -17868, -17699, -17530, -17360, -17189, -17017,
-16845, -16672, -16499, -16325, -16150, -15975, -15799, -15623,
-15446, -15268, -15090, -14911, -14732, -14552, -14372, -14191,
-14009, -13827, -13645, -13462, -13278, -13094, -12909, -12724,
-12539, -12353, -12166, -11980, -11792, -11604, -11416, -11227,
-11038, -10849, -10659, -10469, -10278, -10087, -9895, -9703,
-9511, -9319, -9126, -8932, -8739, -8545, -8351, -8156,
-7961, -7766, -7571, -7375, -7179, -6982, -6786, -6589,
-6392, -6195, -5997, -5799, -5601, -5403, -5205, -5006,
-4807, -4608, -4409, -4210, -4011, -3811, -3611, -3411,
-3211, -3011, -2811, -2610, -2410, -2209, -2009, -1808,
-1607, -1406, -1206, -1005, -804, -603, -402, -201,
};
/*
fix_fft() - perform fast Fourier transform.
fr[n],fi[n] are real,imaginary arrays, INPUT AND RESULT.
size of data = 2**m
set inverse to 0=dft, 1=idft
*/
int fix_fft(int[] fr, int[] fi, int m, boolean inverse)
{
int mr,nn,i,j,l,k,istep, n, scale;
int qr,qi,tr,ti,wr,wi,t;
boolean shift;
n = 1< 1024 /*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 = 10 /*LOG2_N_WAVE*/ -1;
while(l < n) {
if(inverse) {
/* variable scaling, depending upon data */
shift = false;
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 = true;
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 = true;
}
/* 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 = Sinewave[j+1024 /*N_WAVE*/ /4];
wi = -Sinewave[j];
if(inverse)
wi = -wi;
if(shift) {
wr >>= 1;
wi >>= 1;
}
for(i=m; i< n; i+=istep) {
j = i + l;
tr = ((wr*fr[j])>>15) - ((wi*fi[j])>>15);
ti = ((wr*fi[j])>>15) + ((wi*fr[j])>>15);
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;
}