diff options
Diffstat (limited to 'soundfft.c')
| -rw-r--r-- | soundfft.c | 196 |
1 files changed, 196 insertions, 0 deletions
diff --git a/soundfft.c b/soundfft.c new file mode 100644 index 0000000..b1a7f81 --- /dev/null +++ b/soundfft.c @@ -0,0 +1,196 @@ +/* + * Program: REALFFT.C + * Author: Philip VanBaren + * Date: 2 September 1993 + * + * Description: These routines perform an FFT on real data. + * On a 486/33 compiled using Borland C++ 3.1 with full + * speed optimization and a small memory model, a 1024 point + * FFT takes about 16ms. + * This code is for integer data, but could be converted + * to float or double simply by changing the data types + * and getting rid of the bit-shifting necessary to prevent + * overflow/underflow in fixed-point calculations. + * + * Note: Output is BIT-REVERSED! so you must use the BitReversed to + * get legible output, (i.e. Real_i = buffer[ BitReversed[i] ] + * Imag_i = buffer[ BitReversed[i]+1 ] ) + * Input is in normal order. + * + * Copyright (C) 1995 Philip VanBaren + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2 of the License, or + * (at your option) any later version. + * + * This program 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 General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. + */ + +/* + * For Borland C++/TASM you may define the flag ASMFFTCODE in order to use the + * assembly code for the FFT. This will skip the compilation of the C version + * of the routine; you must be sure to add realffta.asm to your makefile + */ + +#include <math.h> +#include <stdlib.h> +#include <stdio.h> + +#include "soundfft.h" + +int *BitReversed; +short *SinTable; +int Points = 0; + +/* + * Initialize the Sine table and Twiddle pointers (bit-reversed pointers) + * for the FFT routine. + */ +void InitializeFFT(int fftlen) +{ + int i; + int temp; + int mask; + + /* + * FFT size is only half the number of data points + * The full FFT output can be reconstructed from this FFT's output. + * (This optimization can be made since the data is real.) + */ + Points = fftlen; + + if((SinTable=(short *)malloc(Points*sizeof(short)))==NULL) + { + puts("Error allocating memory for Sine table."); + exit(1); + } + if((BitReversed=(int *)malloc(Points/2*sizeof(int)))==NULL) + { + puts("Error allocating memory for BitReversed."); + exit(1); + } + + for(i=0;i<Points/2;i++) + { + temp=0; + for(mask=Points/4;mask>0;mask >>= 1) + temp=(temp >> 1) + (i&mask ? Points/2 : 0); + + BitReversed[i]=temp; + } + + for(i=0;i<Points/2;i++) + { + register double s,c; + s=floor(-32768.0*sin(2*M_PI*i/(Points))+0.5); + c=floor(-32768.0*cos(2*M_PI*i/(Points))+0.5); + if(s>32767.5) s=32767; + if(c>32767.5) c=32767; + SinTable[BitReversed[i] ]=(short)s; + SinTable[BitReversed[i]+1]=(short)c; + } +} + +/* + * Free up the memory allotted for Sin table and Twiddle Pointers + */ +void EndFFT(void) +{ + free(BitReversed); + free(SinTable); + Points=0; +} + +/* Include this only if you don't use the REALFFTA.ASM routine */ +/* This is for DOS only */ + +#ifndef ASMFFTCODE + +short *A,*B; +short *sptr; +short *endptr1,*endptr2; +int *br1,*br2; +long HRplus,HRminus,HIplus,HIminus; + +/* + * Actual FFT routine. Must call InitializeFFT(fftlen) first! + */ +void RealFFT(short *buffer) +{ + int ButterfliesPerGroup=Points/4; + + endptr1=buffer+Points; + + /* + * Butterfly: + * Ain-----Aout + * \ / + * / \ + * Bin-----Bout + */ + + while(ButterfliesPerGroup>0) + { + A=buffer; + B=buffer+ButterfliesPerGroup*2; + sptr=SinTable; + + while(A<endptr1) + { + register short sin=*sptr; + register short cos=*(sptr+1); + endptr2=B; + while(A<endptr2) + { + long v1=((long)*B*cos + (long)*(B+1)*sin) >> 15; + long v2=((long)*B*sin - (long)*(B+1)*cos) >> 15; + *B=(*A+v1)>>1; + *(A++)=*(B++)-v1; + *B=(*A-v2)>>1; + *(A++)=*(B++)+v2; + } + A=B; + B+=ButterfliesPerGroup*2; + sptr+=2; + } + ButterfliesPerGroup >>= 1; + } + /* + * Massage output to get the output for a real input sequence. + */ + br1=BitReversed+1; + br2=BitReversed+Points/2-1; + + while(br1<=br2) + { + register long temp1,temp2; + short sin=SinTable[*br1]; + short cos=SinTable[*br1+1]; + A=buffer+*br1; + B=buffer+*br2; + HRplus = (HRminus = *A - *B ) + (*B << 1); + HIplus = (HIminus = *(A+1) - *(B+1)) + (*(B+1) << 1); + temp1 = ((long)sin*HRminus - (long)cos*HIplus) >> 15; + temp2 = ((long)cos*HRminus + (long)sin*HIplus) >> 15; + *B = (*A = (HRplus + temp1) >> 1) - temp1; + *(B+1) = (*(A+1) = (HIminus + temp2) >> 1) - HIminus; + + br1++; + br2--; + } + /* + * Handle DC bin separately + */ + buffer[0]+=buffer[1]; + buffer[1]=0; +} +#endif + |