Commit af810089 authored by Gargaj's avatar Gargaj
Browse files

Merge branch 'drop-the-bass'

parents 77eb00a9 b1825640
#!/usr/bin/env python
# Copyright (c) 2003-2010, Mark Borgerding. All rights reserved.
# This file is part of KISS FFT - https://github.com/mborgerding/kissfft
#
# SPDX-License-Identifier: BSD-3-Clause
# See COPYING file for more information.
import math
import sys
import os
import random
import struct
import popen2
import getopt
import numpy
pi=math.pi
e=math.e
j=complex(0,1)
doreal=0
datatype = os.environ.get('DATATYPE','float')
util = '../tools/fft_' + datatype
minsnr=90
if datatype == 'double':
fmt='d'
elif datatype=='int16_t':
fmt='h'
minsnr=10
elif datatype=='int32_t':
fmt='i'
elif datatype=='simd':
fmt='4f'
sys.stderr.write('testkiss.py does not yet test simd')
sys.exit(0)
elif datatype=='float':
fmt='f'
else:
sys.stderr.write('unrecognized datatype %s\n' % datatype)
sys.exit(1)
def dopack(x,cpx=1):
x = numpy.reshape( x, ( numpy.size(x),) )
if cpx:
s = ''.join( [ struct.pack(fmt*2,c.real,c.imag) for c in x ] )
else:
s = ''.join( [ struct.pack(fmt,c.real) for c in x ] )
return s
def dounpack(x,cpx):
uf = fmt * ( len(x) / struct.calcsize(fmt) )
s = struct.unpack(uf,x)
if cpx:
return numpy.array(s[::2]) + numpy.array( s[1::2] )*j
else:
return numpy.array(s )
def make_random(dims=[1]):
res = []
for i in range(dims[0]):
if len(dims)==1:
r=random.uniform(-1,1)
if doreal:
res.append( r )
else:
i=random.uniform(-1,1)
res.append( complex(r,i) )
else:
res.append( make_random( dims[1:] ) )
return numpy.array(res)
def flatten(x):
ntotal = numpy.size(x)
return numpy.reshape(x,(ntotal,))
def randmat( ndims ):
dims=[]
for i in range( ndims ):
curdim = int( random.uniform(2,5) )
if doreal and i==(ndims-1):
curdim = int(curdim/2)*2 # force even last dimension if real
dims.append( curdim )
return make_random(dims )
def test_fft(ndims):
x=randmat( ndims )
if doreal:
xver = numpy.fft.rfftn(x)
else:
xver = numpy.fft.fftn(x)
open('/tmp/fftexp.dat','w').write(dopack( flatten(xver) , True ) )
x2=dofft(x,doreal)
err = xver - x2
errf = flatten(err)
xverf = flatten(xver)
errpow = numpy.vdot(errf,errf)+1e-10
sigpow = numpy.vdot(xverf,xverf)+1e-10
snr = 10*math.log10(abs(sigpow/errpow) )
print 'SNR (compared to NumPy) : %.1fdB' % float(snr)
if snr<minsnr:
print 'xver=',xver
print 'x2=',x2
print 'err',err
sys.exit(1)
def dofft(x,isreal):
dims=list( numpy.shape(x) )
x = flatten(x)
scale=1
if datatype=='int16_t':
x = 32767 * x
scale = len(x) / 32767.0
elif datatype=='int32_t':
x = 2147483647.0 * x
scale = len(x) / 2147483647.0
cmd='%s -n ' % util
cmd += ','.join([str(d) for d in dims])
if doreal:
cmd += ' -R '
print cmd
p = popen2.Popen3(cmd )
open('/tmp/fftin.dat','w').write(dopack( x , isreal==False ) )
p.tochild.write( dopack( x , isreal==False ) )
p.tochild.close()
res = dounpack( p.fromchild.read() , 1 )
open('/tmp/fftout.dat','w').write(dopack( flatten(res) , True ) )
if doreal:
dims[-1] = int( dims[-1]/2 ) + 1
res = scale * res
p.wait()
return numpy.reshape(res,dims)
def main():
opts,args = getopt.getopt(sys.argv[1:],'r')
opts=dict(opts)
global doreal
doreal = opts.has_key('-r')
if doreal:
print 'Testing multi-dimensional real FFTs'
else:
print 'Testing multi-dimensional FFTs'
for dim in range(1,4):
test_fft( dim )
if __name__ == "__main__":
main()
/*
* Copyright (c) 2003-2010, Mark Borgerding. All rights reserved.
* This file is part of KISS FFT - https://github.com/mborgerding/kissfft
*
* SPDX-License-Identifier: BSD-3-Clause
* See COPYING file for more information.
*/
#include <stdlib.h>
#include <string.h>
#include <stdio.h>
#include "kiss_fft.h"
#include "kiss_fftr.h"
#include <limits.h>
static
double two_tone_test( int nfft, int bin1,int bin2)
{
kiss_fftr_cfg cfg = NULL;
kiss_fft_cpx *kout = NULL;
kiss_fft_scalar *tbuf = NULL;
int i;
double f1 = bin1*2*M_PI/nfft;
double f2 = bin2*2*M_PI/nfft;
double sigpow=0;
double noisepow=0;
#if FIXED_POINT==32
long maxrange = LONG_MAX;
#else
long maxrange = SHRT_MAX;/* works fine for float too*/
#endif
cfg = kiss_fftr_alloc(nfft , 0, NULL, NULL);
tbuf = KISS_FFT_MALLOC(nfft * sizeof(kiss_fft_scalar));
kout = KISS_FFT_MALLOC(nfft * sizeof(kiss_fft_cpx));
/* generate a signal with two tones*/
for (i = 0; i < nfft; i++) {
#ifdef USE_SIMD
tbuf[i] = _mm_set1_ps( (maxrange>>1)*cos(f1*i)
+ (maxrange>>1)*cos(f2*i) );
#else
tbuf[i] = (maxrange>>1)*cos(f1*i)
+ (maxrange>>1)*cos(f2*i);
#endif
}
kiss_fftr(cfg, tbuf, kout);
for (i=0;i < (nfft/2+1);++i) {
#ifdef USE_SIMD
double tmpr = (double)*(float*)&kout[i].r / (double)maxrange;
double tmpi = (double)*(float*)&kout[i].i / (double)maxrange;
#else
double tmpr = (double)kout[i].r / (double)maxrange;
double tmpi = (double)kout[i].i / (double)maxrange;
#endif
double mag2 = tmpr*tmpr + tmpi*tmpi;
if (i!=0 && i!= nfft/2)
mag2 *= 2; /* all bins except DC and Nyquist have symmetric counterparts implied*/
/* if there is power in one of the expected bins, it is signal, otherwise noise*/
if ( i!=bin1 && i != bin2 )
noisepow += mag2;
else
sigpow += mag2;
}
kiss_fft_cleanup();
/*printf("TEST %d,%d,%d noise @ %fdB\n",nfft,bin1,bin2,10*log10(noisepow/sigpow +1e-30) );*/
return 10*log10(sigpow/(noisepow+1e-50) );
}
int main(int argc,char ** argv)
{
int nfft = 4*2*2*3*5;
if (argc>1) nfft = atoi(argv[1]);
int i,j;
double minsnr = 500;
double maxsnr = -500;
double snr;
for (i=0;i<nfft/2;i+= (nfft>>4)+1) {
for (j=i;j<nfft/2;j+=(nfft>>4)+7) {
snr = two_tone_test(nfft,i,j);
if (snr<minsnr) {
minsnr=snr;
}
if (snr>maxsnr) {
maxsnr=snr;
}
}
}
snr = two_tone_test(nfft,nfft/2,nfft/2);
if (snr<minsnr) minsnr=snr;
if (snr>maxsnr) maxsnr=snr;
printf("TwoToneTest: snr ranges from %ddB to %ddB\n",(int)minsnr,(int)maxsnr);
printf("sizeof(kiss_fft_scalar) = %d\n",(int)sizeof(kiss_fft_scalar) );
return 0;
}
WARNINGS=-W -Wall -Wstrict-prototypes -Wmissing-prototypes -Waggregate-return \
-Wcast-align -Wcast-qual -Wnested-externs -Wshadow -Wbad-function-cast \
-Wwrite-strings
ifeq "$(DATATYPE)" ""
DATATYPE=float
endif
ifeq "$(DATATYPE)" "int32_t"
TYPEFLAGS=-DFIXED_POINT=32
endif
ifeq "$(DATATYPE)" "int16_t"
TYPEFLAGS=-DFIXED_POINT=16
endif
ifeq "$(DATATYPE)" "simd"
TYPEFLAGS=-DUSE_SIMD=1 -msse
endif
ifeq "$(TYPEFLAGS)" ""
TYPEFLAGS=-Dkiss_fft_scalar=$(DATATYPE)
endif
ifneq ("$(KISS_FFT_USE_ALLOCA)","")
CFLAGS+= -DKISS_FFT_USE_ALLOCA=1
endif
CFLAGS+= $(CFLAGADD)
FFTUTIL=fft_$(DATATYPE)
FASTFILT=fastconv_$(DATATYPE)
FASTFILTREAL=fastconvr_$(DATATYPE)
PSDPNG=psdpng_$(DATATYPE)
DUMPHDR=dumphdr_$(DATATYPE)
all: $(FFTUTIL) $(FASTFILT) $(FASTFILTREAL)
# $(PSDPNG)
# $(DUMPHDR)
#CFLAGS=-Wall -O3 -pedantic -march=pentiumpro -ffast-math -fomit-frame-pointer $(WARNINGS)
# If the above flags do not work, try the following
CFLAGS=-Wall -O3 $(WARNINGS)
# tip: try -openmp or -fopenmp to use multiple cores
$(FASTFILTREAL): ../kiss_fft.c kiss_fastfir.c kiss_fftr.c
$(CC) -o $@ $(CFLAGS) -I.. $(TYPEFLAGS) -DREAL_FASTFIR $+ -DFAST_FILT_UTIL -lm
$(FASTFILT): ../kiss_fft.c kiss_fastfir.c
$(CC) -o $@ $(CFLAGS) -I.. $(TYPEFLAGS) $+ -DFAST_FILT_UTIL -lm
$(FFTUTIL): ../kiss_fft.c fftutil.c kiss_fftnd.c kiss_fftr.c kiss_fftndr.c
$(CC) -o $@ $(CFLAGS) -I.. $(TYPEFLAGS) $+ -lm
$(PSDPNG): ../kiss_fft.c psdpng.c kiss_fftr.c
$(CC) -o $@ $(CFLAGS) -I.. $(TYPEFLAGS) $+ -lpng -lm
$(DUMPHDR): ../kiss_fft.c dumphdr.c
$(CC) -o $@ $(CFLAGS) -I.. $(TYPEFLAGS) $+ -lm
clean:
rm -f *~ fft fft_* fastconv fastconv_* fastconvr fastconvr_* psdpng psdpng_*
This diff is collapsed.
/*
* Copyright (c) 2003-2004, Mark Borgerding. All rights reserved.
* This file is part of KISS FFT - https://github.com/mborgerding/kissfft
*
* SPDX-License-Identifier: BSD-3-Clause
* See COPYING file for more information.
*/
#include "kfc.h"
typedef struct cached_fft *kfc_cfg;
struct cached_fft
{
int nfft;
int inverse;
kiss_fft_cfg cfg;
kfc_cfg next;
};
static kfc_cfg cache_root=NULL;
static int ncached=0;
static kiss_fft_cfg find_cached_fft(int nfft,int inverse)
{
size_t len;
kfc_cfg cur=cache_root;
kfc_cfg prev=NULL;
while ( cur ) {
if ( cur->nfft == nfft && inverse == cur->inverse )
break;/*found the right node*/
prev = cur;
cur = prev->next;
}
if (cur== NULL) {
/* no cached node found, need to create a new one*/
kiss_fft_alloc(nfft,inverse,0,&len);
#ifdef USE_SIMD
int padding = (16-sizeof(struct cached_fft)) & 15;
// make sure the cfg aligns on a 16 byte boundary
len += padding;
#endif
cur = (kfc_cfg)KISS_FFT_MALLOC((sizeof(struct cached_fft) + len ));
if (cur == NULL)
return NULL;
cur->cfg = (kiss_fft_cfg)(cur+1);
#ifdef USE_SIMD
cur->cfg = (kiss_fft_cfg) ((char*)(cur+1)+padding);
#endif
kiss_fft_alloc(nfft,inverse,cur->cfg,&len);
cur->nfft=nfft;
cur->inverse=inverse;
cur->next = NULL;
if ( prev )
prev->next = cur;
else
cache_root = cur;
++ncached;
}
return cur->cfg;
}
void kfc_cleanup(void)
{
kfc_cfg cur=cache_root;
kfc_cfg next=NULL;
while (cur){
next = cur->next;
free(cur);
cur=next;
}
ncached=0;
cache_root = NULL;
}
void kfc_fft(int nfft, const kiss_fft_cpx * fin,kiss_fft_cpx * fout)
{
kiss_fft( find_cached_fft(nfft,0),fin,fout );
}
void kfc_ifft(int nfft, const kiss_fft_cpx * fin,kiss_fft_cpx * fout)
{
kiss_fft( find_cached_fft(nfft,1),fin,fout );
}
#ifdef KFC_TEST
static void check(int nc)
{
if (ncached != nc) {
fprintf(stderr,"ncached should be %d,but it is %d\n",nc,ncached);
exit(1);
}
}
int main(void)
{
kiss_fft_cpx buf1[1024],buf2[1024];
memset(buf1,0,sizeof(buf1));
check(0);
kfc_fft(512,buf1,buf2);
check(1);
kfc_fft(512,buf1,buf2);
check(1);
kfc_ifft(512,buf1,buf2);
check(2);
kfc_cleanup();
check(0);
return 0;
}
#endif
/*
* Copyright (c) 2003-2004, Mark Borgerding. All rights reserved.
* This file is part of KISS FFT - https://github.com/mborgerding/kissfft
*
* SPDX-License-Identifier: BSD-3-Clause
* See COPYING file for more information.
*/
#ifndef KFC_H
#define KFC_H
#include "kiss_fft.h"
#ifdef __cplusplus
extern "C" {
#endif
/*
KFC -- Kiss FFT Cache
Not needing to deal with kiss_fft_alloc and a config
object may be handy for a lot of programs.
KFC uses the underlying KISS FFT functions, but caches the config object.
The first time kfc_fft or kfc_ifft for a given FFT size, the cfg
object is created for it. All subsequent calls use the cached
configuration object.
NOTE:
You should probably not use this if your program will be using a lot
of various sizes of FFTs. There is a linear search through the
cached objects. If you are only using one or two FFT sizes, this
will be negligible. Otherwise, you may want to use another method
of managing the cfg objects.
There is no automated cleanup of the cached objects. This could lead
to large memory usage in a program that uses a lot of *DIFFERENT*
sized FFTs. If you want to force all cached cfg objects to be freed,
call kfc_cleanup.
*/
/*forward complex FFT */
void kfc_fft(int nfft, const kiss_fft_cpx * fin,kiss_fft_cpx * fout);
/*reverse complex FFT */
void kfc_ifft(int nfft, const kiss_fft_cpx * fin,kiss_fft_cpx * fout);
/*free all cached objects*/
void kfc_cleanup(void);
#ifdef __cplusplus
}
#endif
#endif
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
......@@ -2,7 +2,9 @@
namespace FFT
{
extern float fAmplification;
bool Open();
bool GetFFT( float * samples );
bool GetFFT( float * _samples );
void Close();
}
\ No newline at end of file
This diff is collapsed.
This diff is collapsed.
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment