#include "kiss_fft.h"
#include "kiss_fft2d.h"
#include "_kiss_fft_guts.h"

#define xstr(s) str(s)
#define str(s) #s
double snr_compare( kiss_fft_cpx * test_vec_out,kiss_fft_cpx * testbuf, int n)
{
    int k;
    double sigpow,noisepow,err,snr,scale=0;
    sigpow = noisepow = .000000000000000000000000000001; 

    for (k=0;k<n;++k) {
        sigpow += test_vec_out[k].r * test_vec_out[k].i + 
                  test_vec_out[k].i * test_vec_out[k].i;
        err = test_vec_out[k].r - testbuf[k].r;
        noisepow += err * err;
        err = test_vec_out[k].i - testbuf[k].i;
        noisepow += err * err;

        if (test_vec_out[k].r)
            scale += testbuf[k].r / test_vec_out[k].r;
    }
    snr = 10*log10( sigpow / noisepow );
    scale /= n;
    if (snr<10)
        printf( "\npoor snr, try a scaling factor %f\n" , scale );
    return snr;
}

int test_stride();
int test2d();

int main() { int exit_code=0;

#define NFFT 10
    {
        double snr;
        kiss_fft_cpx test_vec_in[NFFT] = { {14733 ,-20989 },{26000 ,2495 },{-22357 ,-3847 },{27716 ,12164 },{-31059 ,139 },{19845 ,19383 },{8812 ,-6250 },{-13047 ,21227 },{5645 ,15891 },{16370 ,17397 }};
        kiss_fft_cpx test_vec_out[NFFT] = {{52658 ,57610 },{5071.7942464632433257065713405609130859375 ,-20292.93526405603552120737731456756591796875 },{15925.6188406528017367236316204071044921875 ,-40764.3196791790760471485555171966552734375 },{9140.470103009121885406784713268280029296875 ,5124.34458797621846315450966358184814453125 },{15304.821807576650826376862823963165283203125 ,-93210.235176821399363689124584197998046875 },{-101110 ,-87721.999999999898136593401432037353515625 },{20028.8367977931993664242327213287353515625 ,95780.48305671053822152316570281982421875 },{-8524.447753100639602052979171276092529296875 ,-67398.910839595715515315532684326171875 },{68972.722553977349889464676380157470703125 ,-27445.92820071005917270667850971221923828125 },{69862.183403628281666897237300872802734375 ,-31570.49848432456565205939114093780517578125 }};
        kiss_fft_cpx testbuf[NFFT];
        void * cfg = kiss_fft_alloc(NFFT,0,0,0);

        kiss_fft(cfg,test_vec_in,testbuf);
        snr = snr_compare(test_vec_out,testbuf,NFFT);
        printf("DATATYPE=" xstr(kiss_fft_scalar) ", FFT n=%d, inverse=%d, snr = %g dB\n",NFFT,0,snr);
        if (snr<100)
            exit_code++;
        free(cfg);
    }
#undef NFFT    

#define NFFT 10
    {
        double snr;
        kiss_fft_cpx test_vec_in[NFFT] = { {2763 ,8156 },{15120 ,-22117 },{23990 ,22942 },{19376 ,27663 },{-28261 ,-19449 },{10065 ,-143 },{19443 ,-15224 },{24506 ,-30762 },{-10631 ,-30800 },{-8842 ,1465 }};
        kiss_fft_cpx test_vec_out[NFFT] = {{67529 ,-58269 },{-94853.81692500403732992708683013916015625 ,32264.0256640331936068832874298095703125 },{-13103.645618317759726778604090213775634765625 ,91304.50185425896779634058475494384765625 },{105102.842114731451147235929965972900390625 ,-32097.974087137539754621684551239013671875 },{39495.306403543901978991925716400146484375 ,53704.736317971357493661344051361083984375 },{-52920.999999999898136593401432037353515625 ,-10481.00000000006184563972055912017822265625 },{25647.3834779135140706785023212432861328125 ,45067.195694883936084806919097900390625 },{-79649.431770286828395910561084747314453125 ,47730.743283336036256514489650726318359375 },{-55428.0442631396581418812274932861328125 ,-91742.43386711427592672407627105712890625 },{85811.40658055929816327989101409912109375 ,4079.205139768380831810645759105682373046875 }};
        kiss_fft_cpx testbuf[NFFT];
        void * cfg = kiss_fft_alloc(NFFT,1,0,0);

        kiss_fft(cfg,test_vec_in,testbuf);
        snr = snr_compare(test_vec_out,testbuf,NFFT);
        printf("DATATYPE=" xstr(kiss_fft_scalar) ", FFT n=%d, inverse=%d, snr = %g dB\n",NFFT,1,snr);
        if (snr<100)
            exit_code++;
        free(cfg);
    }
#undef NFFT    

#define NFFT 12
    {
        double snr;
        kiss_fft_cpx test_vec_in[NFFT] = { {3044 ,-10605 },{24938 ,-27647 },{23216 ,30804 },{6706 ,-3134 },{-14653 ,28784 },{-18566 ,-14839 },{-9220 ,-28786 },{-9510 ,-5485 },{4530 ,8300 },{24309 ,-22830 },{23926 ,24160 },{-18558 ,-1011 }};
        kiss_fft_cpx test_vec_out[NFFT] = {{40162 ,-22289 },{95930.808466660906560719013214111328125 ,37514.387642344343475997447967529296875 },{-93500.94261657944298349320888519287109375 ,-145451.45430378968012519180774688720703125 },{-100687 ,-51742 },{14196.6405601739315898157656192779541015625 ,-116401.279735569914919324219226837158203125 },{-10735.80846666081924922764301300048828125 ,17493.61235765566016198135912418365478515625 },{21524 ,127603 },{32849.26584106768132187426090240478515625 ,51183.299072623354732058942317962646484375 },{20158.3594398259956506080925464630126953125 ,-57374.720264430143288336694240570068359375 },{10685 ,52344 },{-39596.057383420717087574303150177001953125 ,-22432.545696210218011401593685150146484375 },{45541.7341589324714732356369495391845703125 ,2292.70092737659797421656548976898193359375 }};
        kiss_fft_cpx testbuf[NFFT];
        void * cfg = kiss_fft_alloc(NFFT,0,0,0);

        kiss_fft(cfg,test_vec_in,testbuf);
        snr = snr_compare(test_vec_out,testbuf,NFFT);
        printf("DATATYPE=" xstr(kiss_fft_scalar) ", FFT n=%d, inverse=%d, snr = %g dB\n",NFFT,0,snr);
        if (snr<100)
            exit_code++;
        free(cfg);
    }
#undef NFFT    

#define NFFT 12
    {
        double snr;
        kiss_fft_cpx test_vec_in[NFFT] = { {29648 ,4570 },{13259 ,-2954 },{19399 ,-17312 },{-25117 ,18529 },{-32521 ,-26069 },{-4116 ,3672 },{30632 ,32625 },{-20086 ,28211 },{-19595 ,2048 },{-1288 ,16246 },{-3780 ,-22556 },{20756 ,12759 }};
        kiss_fft_cpx test_vec_out[NFFT] = {{7191 ,49769 },{120952.49999288018443621695041656494140625 ,-65811.75039975097752176225185394287109375 },{73304.9862905458430759608745574951171875 ,66153.32905302778817713260650634765625 },{-26184 ,20094 },{68464.933281851219362579286098480224609375 ,31479.63464575339457951486110687255859375 },{-19499.4999928800607449375092983245849609375 ,-45330.2496002491025137715041637420654296875 },{40375 ,-103157 },{-15568.3018713572455453686416149139404296875 ,11613.467329754619640880264341831207275390625 },{25969.0667181486569461412727832794189453125 ,134661.3653542466345243155956268310546875 },{-111254 ,-44510 },{146375.01370945409871637821197509765625 ,44263.670946972328238189220428466796875 },{45649.3018713572746491990983486175537109375 ,-44385.4673297546760295517742633819580078125 }};
        kiss_fft_cpx testbuf[NFFT];
        void * cfg = kiss_fft_alloc(NFFT,1,0,0);

        kiss_fft(cfg,test_vec_in,testbuf);
        snr = snr_compare(test_vec_out,testbuf,NFFT);
        printf("DATATYPE=" xstr(kiss_fft_scalar) ", FFT n=%d, inverse=%d, snr = %g dB\n",NFFT,1,snr);
        if (snr<100)
            exit_code++;
        free(cfg);
    }
#undef NFFT    

#define NFFT 14
    {
        double snr;
        kiss_fft_cpx test_vec_in[NFFT] = { {-17717 ,-15600 },{-17788 ,19203 },{22838 ,13954 },{26984 ,-28478 },{-3626 ,-29644 },{6324 ,-23085 },{15432 ,-22201 },{-10999 ,800 },{-3576 ,25838 },{30176 ,-7785 },{-1584 ,-21386 },{24920 ,-2334 },{-18853 ,-8579 },{15926 ,14105 }};
        kiss_fft_cpx test_vec_out[NFFT] = {{68457 ,-85192 },{-73157.08577897600480355322360992431640625 ,29842.11742827893976937048137187957763671875 },{-2895.2606216588173992931842803955078125 ,63183.58866548654623329639434814453125 },{-27915.0677076224164920859038829803466796875 ,-48444.4766526431267266161739826202392578125 },{11071.092839997043483890593051910400390625 ,31969.9043844794505275785923004150390625 },{4977.4895202731959216180257499217987060546875 ,4130.75877337983911274932324886322021484375 },{-46830.001595729423570446670055389404296875 ,24345.1362955485892598517239093780517578125 },{-82628.999999999927240423858165740966796875 ,-30043.9999999999927240423858165740966796875 },{1092.62014028533667442388832569122314453125 ,-115886.35438792980858124792575836181640625 },{116078.984861459306557662785053253173828125 ,-53628.8559552718434133566915988922119140625 },{-87675.824885991038172505795955657958984375 ,-134126.87834761265548877418041229248046875 },{-4256.0042434279157532728277146816253662109375 ,-61673.352946727405651472508907318115234375 },{-144231.62587690286454744637012481689453125 ,112106.60339002744876779615879058837890625 },{19874.6833482937872759066522121429443359375 ,45017.8093529834950459189713001251220703125 }};
        kiss_fft_cpx testbuf[NFFT];
        void * cfg = kiss_fft_alloc(NFFT,0,0,0);

        kiss_fft(cfg,test_vec_in,testbuf);
        snr = snr_compare(test_vec_out,testbuf,NFFT);
        printf("DATATYPE=" xstr(kiss_fft_scalar) ", FFT n=%d, inverse=%d, snr = %g dB\n",NFFT,0,snr);
        if (snr<100)
            exit_code++;
        free(cfg);
    }
#undef NFFT    

#define NFFT 14
    {
        double snr;
        kiss_fft_cpx test_vec_in[NFFT] = { {18387 ,17443 },{-24633 ,28445 },{-14413 ,-18313 },{19595 ,27925 },{-16406 ,-1049 },{9002 ,-12125 },{-16679 ,-9081 },{-16173 ,-14355 },{2223 ,-3105 },{24907 ,2448 },{4568 ,27578 },{3679 ,-820 },{11877 ,20660 },{-47 ,-21552 }};
        kiss_fft_cpx test_vec_out[NFFT] = {{5887 ,44099 },{33087.4744857999103260226547718048095703125 ,-176.7023630800031241960823535919189453125 },{-84337.5165114859701134264469146728515625 ,-45466.202714726721751503646373748779296875 },{23636.47562663447388331405818462371826171875 ,-31899.788673282004310749471187591552734375 },{-30514.45509819767539738677442073822021484375 ,14432.78612769579558516852557659149169921875 },{-5992.53789459774634451605379581451416015625 ,22688.016220844045164994895458221435546875 },{-44805.7131123236395069397985935211181640625 ,32970.175651043697143904864788055419921875 },{-26772.99999999979627318680286407470703125 ,24166.99999999986539478413760662078857421875 },{153786.34064676126581616699695587158203125 ,-50266.255254303789115510880947113037109375 },{167350.84601866011507809162139892578125 ,12803.59263587257373728789389133453369140625 },{10087.87980353427337831817567348480224609375 ,74187.52814839812344871461391448974609375 },{62821.089251506768050603568553924560546875 ,81601.98288412831607274711132049560546875 },{5394.4642717120805173180997371673583984375 ,-48341.031958107007085345685482025146484375 },{-12210.34748800349552766419947147369384765625 ,113401.899295517083373852074146270751953125 }};
        kiss_fft_cpx testbuf[NFFT];
        void * cfg = kiss_fft_alloc(NFFT,1,0,0);

        kiss_fft(cfg,test_vec_in,testbuf);
        snr = snr_compare(test_vec_out,testbuf,NFFT);
        printf("DATATYPE=" xstr(kiss_fft_scalar) ", FFT n=%d, inverse=%d, snr = %g dB\n",NFFT,1,snr);
        if (snr<100)
            exit_code++;
        free(cfg);
    }
#undef NFFT    

    if ( test_stride() ) 
        exit_code++;
    if ( test2d() ) 
        exit_code++;

    return exit_code;
}

int test_stride()
{
#define SKIP_FACTOR 7
#define FFT_SIZE 1800
    void *cfg;
    kiss_fft_cpx buf1in[FFT_SIZE],buf1out[FFT_SIZE];
    kiss_fft_cpx buf2in[SKIP_FACTOR*FFT_SIZE],buf2out[FFT_SIZE];
    int i;
    memset(buf2in,0,sizeof(buf2in));
    for (i=0;i<FFT_SIZE;++i) {
        buf1in[i].r = rand();
        buf1in[i].i = rand();
        buf2in[SKIP_FACTOR*i] = buf1in[i];
    }
    cfg= kiss_fft_alloc(FFT_SIZE,0,0,0);

    kiss_fft(cfg,buf1in,buf1out);
    kiss_fft_stride(cfg,buf2in,buf2out,SKIP_FACTOR);
    if ( memcmp(buf1out,buf2out,sizeof(buf1out) ) ){
        printf("kiss_fft_stride not working for stride =%d\n",SKIP_FACTOR);
        for (i=0;i<FFT_SIZE;++i) {
            printf(stderr,"good[%d]=",i);pcpx(buf1out+i);
            printf(stderr,"bad [%d]=",i);pcpx(buf2out+i);
        }
        free(cfg);
        return 1;
    }
    printf("kiss_fft_stride seems to be working (stride = %d)\n",SKIP_FACTOR);
    free(cfg);
    return 0;
}

int test2d()
{
    void *cfg;
    kiss_fft_cpx in[3*9]={
        {1,0},{0,0},{0,0},{0,0},{0,0},{0,0},{3,0},{0,0},{0,0},
        {0,0},{0,0},{0,0},{0,0},{5,0},{0,0},{0,0},{0,0},{0,0},
        {0,0},{0,0},{0,0},{0,0},{0,0},{0,0},{0,0},{0,0},{7,0} };
    kiss_fft_cpx out[3*9];
    kiss_fft_cpx ver[3*9]= {
    {16.00000,0.00000},{0.16385,5.38749},{4.54576,7.50952},{-2.00000,1.73205},{-6.20961,9.91626},{-6.20961,-9.91626},{-2.00000,-1.73205},{4.54576,-7.50952},{0.16385,-5.38749},{-2.00000,1.73205},{-6.20961,9.91626},{-6.20961,-9.91626},{-2.00000,-1.73205},{4.54576,-7.50952},{0.16385,-5.38749},{16.00000,0.00000},{0.16385,5.38749},{4.54576,7.50952},{-2.00000,-1.73205},{4.54576,-7.50952},{0.16385,-5.38749},{16.00000,-0.00000},{0.16385,5.38749},{4.54576,7.50952},{-2.00000,1.73205},{-6.20961,9.91626},{-6.20961,-9.91626} };
    
    cfg = kiss_fft2d_alloc(3,9,0,0,0);
    kiss_fft2d(cfg,in,out);
    free(cfg);

    if ( snr_compare( out,ver, 3*9) < 100 ) {
        return 1;
    }
    printf("2d complex fft seems to be working\n");
    return 0;
}