Post by thread:[PIC] FFT integer routine for PIC (if it's CCS-C

Hi Scott , Claudio & list. First, thanks for answer. Scott , what we need is a 8 bar spectrum analizer, speed is not a real problem (for now) but the problem is that floating point requires, let's say 16 real & 16 imaginary array, each of 32 bit, this and the cosine calculation, kills the ram memory even using a 16F876. Claudio, you are right, that was my though too, the problem is the cosine that CCS does it by calculation and eats about 60 memory positions, that and the array is too much for the poor pic :) so I found a routine that does 16 bit Integer FFT and uses a Cosine table, this routine works even great that the floating point one, excelent results in Borland C, but porting to pic present 2 problems. It requires 16 x 16 int mutliply (intermediate 32 bit result), when the CCS cannot handle 32 bit integer (at least the one i have, doesn't), but it was solved doing some shifting, and worked great, no problem with that right now. The other problem (when I'm stuck right now) was the Cosine table, it 's too big, 1024 points, and due to the indirect addressing limitations of the pic (8 bit pointer) i can't use more than 256 element table, also the table has a strange, and exagerate amount of data. I tried to cut the data to 256 element's but don't work, I'm not a math guru at all, and the inisdes of the fft are still a mystery for me :), however I'm beginning to understand a bit, but not enough to solve the cosine table by myself :(( one of the things that i don't understand why is the line "Sinewave[j+N_WAVE/4]" testing another sine/cosine table (below), replaced this by i_sin(j+NWAVE/4) and the result is any kind of crap. so tried replacing for i_sin(j) and the one with the minus ahead replaced by i_cos(j) and it works partially but not really good, have a lot of random noise in the background with this approach. I'm lost, that's why i decide to look for another fft in the mean time or perhaps another cosine table that gives me no such errors. but without the fully understanding of the function itself it's very difficult to me to simply figure it out by myself. can you give me some suggestion?, approach?, a way to... :) or any ideas?. thanks for your time. see you Leandro _____________________cosine routine & integer FFT below____ Here is the cosine table that I was trying to use instead of the 1024 one, integrated in the routine (see more below). ------------------------------------------------------------- int i_sin(int x) { static long int const i_sintable[33] = { 0,1608,3212,4808,6393,7962,9512,11039, 12540,14010,15447,16846,18205,19520,20788,22006, 23170,24279,25330,26320,27246,28106,28899,29622, 30274,30853,31357,31786,32138,32413,32610,32729, 32767 }; char i,qtr=0; int result; if(x<0) { x=-x; qtr++; } if(((int)x)>=16384) x=(int)32767-x; i=(char)(x>>9); x=(x&0x1ff)<<3; result=i_sintable[i]+(int)(mul((i_sintable[i+1]-i_sintable[i])<<2,x)>>8); if(qtr) result = -result; return result; } int i_cos(int x) { return i_sin(x+16384); } ----------------------------------------------------------------- Here is the FFT that I'm testing: (works excellent in Borland C 4 ) ----------------------------------------------------------------- /* fix_fft.c - Fixed-point Fast Fourier Transform */ #define N 16 #define TWOPI 6.283185307179586476925287 #define FS 3900 #define W 2 extern int real[N],imag[N]; /* FIX_MPY() - fixed-point multiplication macro. This macro is a statement, not an expression (uses asm). BEWARE: make sure _DX is not clobbered by evaluating (A) or DEST. args are all of type fixed. Scaling ensures that 32767*32767 = 32767. */ #define FIX_MPY(DEST,A,B) DEST = ((long)(A) * (long)(B))>>15 #define N_WAVE 1024 /* dimension of Sinewave[] */ #define LOG2_N_WAVE 10 /* log2(N_WAVE) */ #ifndef fixed #define fixed int #endif extern fixed Sinewave[N_WAVE]; /* placed at end of this file for clarity */ fixed fix_mpy(fixed a, fixed b); /* fix_fft() - perform fast Fourier transform. if n>0 FFT is done, if n<0 inverse FFT is done 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 m, int inverse) { int mr,nn,i,j,l,k,istep, n, scale, shift; fixed qr,qi,tr,ti,wr,wi,t; n = 1<<m; 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 = real[m]; real[m] = real[mr]; real[mr] = tr; ti = imag[m]; imag[m] = imag[mr]; imag[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 = real[i]; if(j < 0) j = -j; m = imag[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 = Sinewave[j+N_WAVE/4]; wi = -Sinewave[j]; //wr = -i_sin(j); //wi = i_cos(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,real[j]) - fix_mpy(wi,imag[j]); ti = fix_mpy(wr,imag[j]) + fix_mpy(wi,real[j]); //tr = ((wr * real[j])>>15) - ((wi * imag[j])>>15); //ti = ((wr * imag[j])>>15) + ((wi * real[j])>>15); qr = real[i]; qi = imag[i]; if(shift) { qr >>= 1; qi >>= 1; } real[j] = qr - tr; imag[j] = qi - ti; real[i] = qr + tr; imag[i] = qi + ti; } } --k; l = istep; } return scale; } /* window() - apply a Hanning window */ void wind(void) { int i,j,k; int n=4; j = N_WAVE/n; n >>= 1; for(i=0,k=N_WAVE/4; i<n; ++i,k+=j) //real[i] = ((real[i]>>7) * ((16384-(Sinewave[k]>>1)>>7))>>8); FIX_MPY(real[i],real[i],16384-(Sinewave[k]>>1)); n <<= 1; for(k-=j; i<n; ++i,k-=j) //real[i] = ((real[i]>>7) * ((16384-(Sinewave[k]>>1)>>7))>>8); FIX_MPY(real[i],real[i],16384-(Sinewave[k]>>1)); } /* fix_mpy() - fixed-point multiplication */ fixed fix_mpy(fixed a, fixed b) { int c,d,e,f; d=a; e=b; a= ((float)d/256) * ((float)e/128); //printf ("\n\rA= %d B= %d Resultado= %d\n\r",d,e,c); //c= (long)(a) * (long)(b)>>15; //printf ("\n\rA= %d B= %d Resultado= %d\n\r",a,b,c); //exit(0); return a; } #if N_WAVE != 1024 ERROR: N_WAVE != 1024 #endif fixed Sinewave[1024] = { 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, }; -- http://www.piclist.com hint: The PICList is archived three different ways. See http://www.piclist.com/#archives for details.

On Wed, 7 Aug 2002, Leandro wrote: > Hi Scott , Claudio & list. > > First, thanks for answer. > > Scott , what we need is a 8 bar spectrum analizer, speed is not a real > problem (for > now) but the problem is that floating point requires, let's say 16 real & > 16 imaginary > array, each of 32 bit, this and the cosine calculation, kills the ram > memory even using > a 16F876. Okay. {Quote hidden}> > Claudio, you are right, that was my though too, the problem is the cosine > that CCS does it by calculation and eats about 60 memory positions, > that and the array is too much for the poor pic :) > so I found a routine that does 16 bit Integer FFT and uses a Cosine table, > this routine works even great that the floating point one, excelent results > in Borland C, but porting to pic present 2 problems. > > It requires 16 x 16 int mutliply (intermediate 32 bit result), when the CCS > cannot handle 32 bit integer (at least the one i have, doesn't), but it was > solved doing some shifting, and worked great, no problem with that right now. > > The other problem (when I'm stuck right now) was the Cosine table, it 's > too big, 1024 points, and due to the indirect addressing limitations of > the pic (8 bit pointer) i can't use more than 256 element table, also > the table has a strange, and exagerate amount of data. Take a look at: http://www.dattalo.com/technical/software/pic/picsine.html {Quote hidden}> > I tried to cut the data to 256 element's but don't work, I'm not a math guru > at all, and the inisdes of the fft are still a mystery for me :), however I'm > beginning to understand a bit, but not enough to solve the cosine table > by myself :(( > one of the things that i don't understand why is the line > "Sinewave[j+N_WAVE/4]" > testing another sine/cosine table (below), replaced this > by i_sin(j+NWAVE/4) and the result is any kind of crap. > so tried replacing for i_sin(j) and the one with the minus ahead > replaced by i_cos(j) and it works partially but not really good, have > a lot of random noise in the background with this approach. > > I'm lost, that's why i decide to look for another fft in the mean time > or perhaps another cosine table that gives me no such errors. > but without the fully understanding of the function itself it's very > difficult to me to simply figure it out by myself. > > can you give me some suggestion?, approach?, a way to... :) or > any ideas?. If you need *only* 8 frequency banks then an FFT is really a waste of time in my opinion. The difference between an FFT with 8 banks and 8 discrete filters is minimal. Furthermore, for the PIC architecture, the incurred overhead of an FFT computation will rapdily consume the minimal savings that an FFT will provide. Recall that an FFT is an efficient way of computing the DFT (discrete Fourier Transform). The DFT is the dot product between your signal and the complex exponential exp(j*omega*t). In other words, for one frequency, omega_1, you create a vector of equal-spaced samples of exp(j*omega_1*k/N). This vector, or array if you prefer, is multiplied element-by-element by the samples of your input signal. You may wish to experiment with DFT's nd verify that they produce the same result as FFT's. One other problem with FFT's is that the frequency bins are equally spaced. I doubt that this what you really want for your application. At least with a DFT you can select which 8 frequencies you wish to filter. You may also wish to take a look at: http://www.dattalo.com/technical/theory/dtmf.html Where a theoretical technique is presented for decoding DTMF signals using square waves instead of sine waves. For DTMF, there are 8 frequencies like in your case. However, the midpoints of these frequencies are not harmonically related. Scott -- http://www.piclist.com hint: The PICList is archived three different ways. See http://www.piclist.com/#archives for details.