mirror/kissfft
1
0
Fork 0
mirror of https://github.com/mborgerding/kissfft.git synced 2025-06-04 01:28:23 -04:00
Code Issues Packages Projects Releases Wiki Activity

Initial revision

Browse Source
This commit is contained in:
Mark Borgerding 2003-08-09 00:59:17 +00:00
parent 4a113bc191
commit fa03256dc2
5 changed files with 367 additions and 0 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

11
COPYING Normal file
View File

@ -0,0 +1,11 @@
Copyright (c) 2003, Mark Borgerding
All rights reserved.
Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:
* Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.
* Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution.
* Neither the author nor the names of any contributors may be used to endorse or promote products derived from this software without specific prior written permission.
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

20
Makefile Normal file
View File

@ -0,0 +1,20 @@
all: kiss_fft_s kiss_fft_f kiss_fft_d
kiss_fft_s: kiss_fft.h kiss_fft.c
gcc -Wall -O3 -o kiss_fft_s -DFIXED_POINT -DFFT_UTIL kiss_fft.c -lm
kiss_fft_f: kiss_fft.h kiss_fft.c
gcc -Wall -O3 -o kiss_fft_f -Dkiss_fft_scalar=float -DFFT_UTIL kiss_fft.c -lm
kiss_fft_d: kiss_fft.h kiss_fft.c
gcc -Wall -O3 -o kiss_fft_d -Dkiss_fft_scalar=double -DFFT_UTIL kiss_fft.c -lm
clean:
rm -f kiss_fft_s kiss_fft_f kiss_fft_d *~ fftin.dat fftout.dat
test: clean all
./test.oct
tarball: clean
tar -czf kiss_fft.tar.gz .

79
README Normal file
View File

@ -0,0 +1,79 @@
KISS FFT - A power-of-two Fast Fourier Transform based up on the principle,
"Keep It Simple, Stupid."
There are many great fft libraries already around. Kiss FFT is not trying
to be better than any of them. It only attempts to be a reasonably efficient,
moderately useful FFT that can use fixed or floating data types and can be
incorporated into someone's C program in a few minutes with trivial licensing.
USAGE:
The basic usage is:
void * cfg = kiss_fft_alloc( nfft ,inverse_fft );
while ...
...
kiss_fft( cfg , cx_buf_in_out );
...
free(cfg);
Declarations are in "kiss_fft.h", along with a brief description of the
two functions you'll need to use. Code definitions are in kiss_fft.c, along
with sample usage code.
The code can be easily recompiled to work with 16bit fixed point data,
or various floating point types. The default is float.
BACKGROUND:
I started coding this because I couldn't find a fixed point FFT that didn't
use assembly code. I started with floating point numbers so I could get the
theory straight before working on fixed point issues. In the end, I had a
little bit of code that could be recompiled easily to do ffts with short, float,
or double (other types should be easy too).
Once I got my FFT working, I wanted to get some performance numbers against
a well respected and highly optimized fft library. I don't want to criticize
this great library, so let's call it FFT_BRANDX.
During this process, I learned:
1. FFT_BRANDX has 500 times as many lines of code as Kiss
(and that's just the C code).
2. It took me an embarrassingly long time to get FFT_BRANDX working.
3. FFT_BRANDX is almost 3 times faster than Kiss
It is wonderful that free, highly optimized libraries like FFT_BRANDX exist.
But such libraries carry a huge burden of complexity necessary to extract every
last bit of performance.
Sometimes simpler is better, even if it's not better.
PERFORMANCE:
(on Athlon XP 2100+, with gcc 2.96, optimization O3, float data type)
Kiss performed 1000 1024-pt ffts in 136 ms. This translates to 7.5 Msamples/s.
Just for comparison, it took md5sum 160 ms to process the same amount of data
DO NOT:
... use Kiss if you need the absolute fastest fft in the world
... use Kiss if you need mixed radix FFTs
... ask me to add features that will bloat the code
UNDER THE HOOD:
Kiss uses a complex-only, frequency decimation, radix 2, in-place FFT. Bit reversed
addressing is corrected as the last step in the transform. No scaling is done.
LICENSE:
BSD, see COPYING for details.
TODO:
Make the fixed point scaling and bit shifts more easily configurable.
Document/revisit the input/output fft scaling
See if the fixed point code can be optimized a little without adding complexity.
AUTHOR:
Mark Borgerding
mark@borgerding.net

212
kiss_fft.c Normal file
View File

@ -0,0 +1,212 @@
/*
Copyright (c) 2003, Mark Borgerding
All rights reserved.
Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:
* Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.
* Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution.
* Neither the author nor the names of any contributors may be used to endorse or promote products derived from this software without specific prior written permission.
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
*/
#include <stdlib.h>
#include <math.h>
#include "kiss_fft.h"
/*
* kiss_fft.h
* defines kiss_fft_scalar as either short or a float type
* and defines
* typedef struct {
* kiss_fft_scalar r;
* kiss_fft_scalar i;
* }kiss_fft_cpx;
*/
typedef struct {
int nfft;
kiss_fft_cpx* twiddle;
int * swap_indices;
int nswap;
}kiss_fft_state;
#ifdef FIXED_POINT
# define TWIDDLE_RANGE ( (1<<15) - 1 )
# define C_ADD(x,a,b) \
do{ (x).r = ( ( (a).r+(b).r +1) >> 1 );\
(x).i = ( ( (a).i+(b).i +1) >> 1 ); }while(0)
# define C_SUB(x,a,b) \
do{ (x).r = ( ( (a).r-(b).r +1) >> 1 );\
(x).i = ( ( (a).i-(b).i +1) >> 1 ); }while(0)
# define C_MUL(m,a,b) \
do{ (m).r = ( ( (a).r*(b.r) - (a).i*(b).i) + (1<<14) ) >> 15;\
(m).i = ( ( (a).r*(b.i) + (a).i*(b).r) + (1<<14) ) >> 15;\
}while(0)
#else // floating point
#define C_ADD(x,a,b) \
do{ (x).r = (a).r+(b).r;\
(x).i = (a).i+(b).i;}while(0)
#define C_SUB(x,a,b) \
do{ (x).r = (a).r-(b).r;\
(x).i = (a).i-(b).i;}while(0)
#define C_MUL(m,a,b) \
do{ (m).r = (a).r*(b.r) - (a).i*(b).i;\
(m).i = (a).r*(b.i) + (a).i*(b).r; }while(0)
#endif
// create the twiddle factors
static
void make_twiddles(kiss_fft_cpx * twiddle,int nfft,int inverse_fft)
{
const double pi=3.14159265358979323846264338327;
int n;
for (n=0;n< (nfft>>1);++n) {
#ifdef FIXED_POINT
twiddle[n].r = (kiss_fft_scalar)(TWIDDLE_RANGE*cos(2*pi*n/nfft)+.5);
twiddle[n].i = -(kiss_fft_scalar)(TWIDDLE_RANGE*sin(2*pi*n/nfft)+.5);
#else
twiddle[n].r = cos(2*pi*n/nfft);
twiddle[n].i = -sin(2*pi*n/nfft);
#endif
if (inverse_fft)
twiddle[n].i *= -1;
}
}
// reverse the bits for numbers from 0 to N-1
static
int bit_reverse(int i,int N)
{
int rev=0;
while ( N >>= 1 ) {
rev = rev*2 + (i&1);
i>>=1;
}
return rev;
}
// make a list of index swaps that need to be done for bit-reversed addressing
static
int make_bit_reverse_indices(int *swap_indices ,int nfft)
{
int n,nswap;
nswap=0;
// set up swap indices to unwrap bit-reversed addressing
for ( n=0; n<nfft; ++n ) {
swap_indices[nswap*2] = bit_reverse(n,nfft);
if ( n < swap_indices[nswap*2] ) {
swap_indices[nswap*2+1] = n;
++nswap;
}
}
return nswap;
}
static
void undo_bit_rev( kiss_fft_cpx * f, const int * swap_indices,int nswap)
{
int i,i0,i1;
kiss_fft_cpx tmp;
// unwrap bit-reverse indexing
for ( i=0 ; i<nswap ; ++i ) {
i0= *swap_indices++;
i1= *swap_indices++;
tmp = f[i0];
f[i0] = f[i1];
f[i1] = tmp;
}
}
// the heart of the fft
static
void fft_work(int N,kiss_fft_cpx *f,const kiss_fft_cpx * twid,int twid_step)
{
N >>= 1;
{ // declare automatic variables in here to reduce stack usage
int n;
kiss_fft_cpx csum,cdiff;
for (n=0;n<N;++n) {
C_ADD(csum, f[n] , f[N+n] );
C_SUB(cdiff, f[n] , f[N+n] );
f[n]=csum;
C_MUL(f[N+n], cdiff, twid[n*twid_step] );
}
}
if (N==1) return;
fft_work(N,f,twid,twid_step*2);
fft_work(N,f+N,twid,twid_step*2);
}
void * kiss_fft_alloc(int nfft,int inverse_fft)
{
kiss_fft_state * st=NULL;
// allocate one large buffer to hold the state and the buffers for twiddles and bit-rev indices
int size = sizeof(kiss_fft_state) + (nfft>>1)*sizeof(kiss_fft_cpx) + nfft*sizeof(int);
st = ( kiss_fft_state *)malloc(size);
if (st) {
st->twiddle = (kiss_fft_cpx*)(st+1);
st->swap_indices = (int*)( st->twiddle + (nfft>>1) );
st->nfft = nfft;
make_twiddles(st->twiddle,nfft,inverse_fft);
st->nswap = make_bit_reverse_indices(st->swap_indices,nfft);
}
return st;
}
void kiss_fft(const void * cfg,kiss_fft_cpx *f)
{
const kiss_fft_state * st = cfg;
fft_work(st->nfft, f, st->twiddle , 1 );
undo_bit_rev(f,st->swap_indices,st->nswap);
}
#ifdef FFT_UTIL
#include <stdio.h>
#include <string.h>
int main(int argc,char ** argv)
{
kiss_fft_cpx * buf=NULL;
void *st;
int nfft=1024;
int inverse_fft=0;
fprintf(stderr,"sizeof(kiss_fft_cpx) == %d\n",sizeof(kiss_fft_cpx) );
if (argc>1 && 0==strcmp("-i",argv[1])) {
inverse_fft = 1;
--argc;++argv;
}
if (argc>1) {
nfft = atoi(argv[1]);
}
buf=(kiss_fft_cpx*)malloc( sizeof(kiss_fft_cpx) * nfft );
st = kiss_fft_alloc( nfft ,inverse_fft );
while ( fread( buf , sizeof(kiss_fft_cpx) , nfft , stdin ) > 0 ) {
kiss_fft( st , buf );
fwrite( buf , sizeof(kiss_fft_cpx) , nfft , stdout );
}
free(st);
free(buf);
return 0;
}
#endif

45
kiss_fft.h Normal file
View File

@ -0,0 +1,45 @@
#ifndef KISS_FFT_H
#define KISS_FFT_H
#ifndef kiss_fft_scalar
# ifdef FIXED_POINT
# define kiss_fft_scalar short
# else
# define kiss_fft_scalar float
# endif
#endif
typedef struct {
kiss_fft_scalar r;
kiss_fft_scalar i;
}kiss_fft_cpx;
/*
* fft_alloc
*
* Initialize a radix 2 FFT (or IFFT) algorithm.
*
* The return value from fft_alloc is a cfg buffer used internally
* by the fft routine.
*
* Call free() on it when done using it to avoid memory leaks.
* */
void* kiss_fft_alloc(int nfft,int inverse_fft);
// free() the state when done using it
/*
* kiss_fft
*
* Perform an in-place FFT on a complex input buffer.
*
* the input buffer should be real(f[0]) , imag(f[0]), ... ,real(f[nfft-1]) , imag(f[nfft-1])
* the the output will be real(F[0]) , imag(F[0]), ... ,real(F[nfft-1]) , imag(F[nfft-1])
*
* */
void kiss_fft( const void* cfg , kiss_fft_cpx *f ); // call for each buffer
// when done with the cfg for a given fft size and direction, simply free it
#define kiss_fft_free free
#endif