text:
# KISS FFT [![Build Status](https://travis-ci.com/mborgerding/kissfft.svg?branch=master)](https://travis-ci.com/mborgerding/kissfft)

KISS FFT - A mixed-radix 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 for 1-d complex FFT is:

```c
    #include "kiss_fft.h"
    kiss_fft_cfg cfg = kiss_fft_alloc( nfft ,is_inverse_fft ,0,0 );
    while ...

        ... // put kth sample in cx_in[k].r and cx_in[k].i

        kiss_fft( cfg , cx_in , cx_out );

        ... // transformed. DC is in cx_out[0].r and cx_out[0].i

    kiss_fft_free(cfg);
```
 - **Note**: frequency-domain data is stored from dc up to 2pi.
    so cx_out[0] is the dc bin of the FFT
    and cx_out[nfft/2] is the Nyquist bin (if exists)

Declarations are in "kiss_fft.h", along with a brief description of the
functions you'll need to use.

Code definitions for 1d complex FFTs are in kiss_fft.c.

You can do other cool stuff with the extras you'll find in tools/
> - multi-dimensional FFTs
> - real-optimized FFTs  (returns the positive half-spectrum:
    (nfft/2+1) complex frequency bins)
> - fast convolution FIR filtering (not available for fixed point)
> - spectrum image creation

The core fft and most tools/ code can be compiled to use float, double,
 Q15 short or Q31 samples. 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 was curious about the speed compared to
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 more than 100K lines of code. The core of kiss_fft is about 500 lines (cpx 1-d).
> 2. It took me an embarrassingly long time to get FFT_BRANDX working.
> 3. A simple program using FFT_BRANDX is 522KB. A similar program using kiss_fft is 18KB (without optimizing for size).
> 4. FFT_BRANDX is roughly twice as fast as KISS FFT in default mode.

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.**

## FREQUENTLY ASKED QUESTIONS:
> Q: Can I use kissfft in a project with a ___ license?

> A: Yes.  See LICENSE below.

> Q: Why don't I get the output I expect?

> A: The two most common causes of this are
> 	1) scaling : is there a constant multiplier between what you got and what you want?
> 	2) mixed build environment -- all code must be compiled with same preprocessor
> 	definitions for FIXED_POINT and kiss_fft_scalar

> Q: Will you write/debug my code for me?

> A: Probably not unless you pay me.  I am happy to answer pointed and topical questions, but
> I may refer you to a book, a forum, or some other resource.

## PERFORMANCE
    (on Athlon XP 2100+, with gcc 2.96, float data type)

Kiss performed 10000 1024-pt cpx ffts in .63 s of cpu time.
For comparison, it took md5sum twice as long to process the same amount of data.
Transforming 5 minutes of CD quality audio takes less than a second (nfft=1024).

**DO NOT:**
- use Kiss if you need the Fastest Fourier Transform in the World
- ask me to add features that will bloat the code

## UNDER THE HOOD

Kiss FFT uses a time decimation, mixed-radix, out-of-place FFT. If you give it an input buffer
and output buffer that are the same, a temporary buffer will be created to hold the data.

No static data is used.  The core routines of kiss_fft are thread-safe (but not all of the tools directory).[

No scaling is done for the floating point version (for speed).
Scaling is done both ways for the fixed-point version (for overflow prevention).

Optimized butterflies are used for factors 2,3,4, and 5.

The real (i.e. not complex) optimization code only works for even length ffts.  It does two half-length
FFTs in parallel (packed into real&imag), and then combines them via twiddling.  The result is
nfft/2+1 complex frequency bins from DC to Nyquist.  If you don't know what this means, search the web.

The fast convolution filtering uses the overlap-scrap method, slightly
modified to put the scrap at the tail.

## LICENSE
    Revised BSD License, see COPYING for verbiage.
    Basically, "free to use&change, give credit where due, no guarantees"
    Note this license is compatible with GPL at one end of the spectrum and closed, commercial software at
    the other end.  See http://www.fsf.org/licensing/licenses

## TODO
 - Add real optimization for odd length FFTs
 - Document/revisit the input/output fft scaling
 - Make doc describing the overlap (tail) scrap fast convolution filtering in kiss_fastfir.c
 - Test all the ./tools/ code with fixed point (kiss_fastfir.c doesn't work, maybe others)

## AUTHOR
    Mark Borgerding
    Mark@Borgerding.net


text:
Speed:
    * If you want to use multiple cores, then compile with -openmp or -fopenmp (see your compiler docs).
	Realize that larger FFTs will reap more benefit than smaller FFTs. This generally uses more CPU time, but
	less wall time.

    * experiment with compiler flags
        Special thanks to Oscar Lesta. He suggested some compiler flags
        for gcc that make a big difference. They shave 10-15% off
        execution time on some systems.  Try some combination of:
                -march=pentiumpro
                -ffast-math
                -fomit-frame-pointer

    * If the input data has no imaginary component, use the kiss_fftr code under tools/.
      Real ffts are roughly twice as fast as complex.

    * If you can rearrange your code to do 4 FFTs in parallel and you are on a recent Intel or AMD machine,
    then you might want to experiment with the USE_SIMD code.  See README.simd

Reducing code size:
    * remove some of the butterflies. There are currently butterflies optimized for radices
        2,3,4,5.  It is worth mentioning that you can still use FFT sizes that contain
        other factors, they just won't be quite as fast.  You can decide for yourself
        whether to keep radix 2 or 4.  If you do some work in this area, let me
        know what you find.

    * For platforms where ROM/code space is more plentiful than RAM,
     consider creating a hardcoded kiss_fft_state. In other words, decide which
     FFT size(s) you want and make a structure with the correct factors and twiddles.

    * Frank van der Hulst offered numerous suggestions for smaller code size and correct operation
    on embedded targets.  "I'm happy to help anyone who is trying to implement KISSFFT on a micro"

    Some of these were rolled into the mainline code base:
        - using long casts to promote intermediate results of short*short multiplication
        - delaying allocation of buffers that are sometimes unused.
    In some cases, it may be desirable to limit capability in order to better suit the target:
        - predefining the twiddle tables for the desired fft size.

repo_language: C

file_name: README

mime_type: text/x-c


text:
/*
 *  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 "_kiss_fft_guts.h"
/* The guts header contains all the multiplication and addition macros that are defined for
 fixed or floating point complex numbers.  It also delares the kf_ internal functions.
 */

static void kf_bfly2(
        kiss_fft_cpx * Fout,
        const size_t fstride,
        const kiss_fft_cfg st,
        int m
        )
{
    kiss_fft_cpx * Fout2;
    kiss_fft_cpx * tw1 = st->twiddles;
    kiss_fft_cpx t;
    Fout2 = Fout + m;
    do{
        C_FIXDIV(*Fout,2); C_FIXDIV(*Fout2,2);

        C_MUL (t,  *Fout2 , *tw1);
        tw1 += fstride;
        C_SUB( *Fout2 ,  *Fout , t );
        C_ADDTO( *Fout ,  t );
        ++Fout2;
        ++Fout;
    }while (--m);
}

static void kf_bfly4(
        kiss_fft_cpx * Fout,
        const size_t fstride,
        const kiss_fft_cfg st,
        const size_t m
        )
{
    kiss_fft_cpx *tw1,*tw2,*tw3;
    kiss_fft_cpx scratch[6];
    size_t k=m;
    const size_t m2=2*m;
    const size_t m3=3*m;

    tw3 = tw2 = tw1 = st->twiddles;

    do {
        C_FIXDIV(*Fout,4); C_FIXDIV(Fout[m],4); C_FIXDIV(Fout[m2],4); C_FIXDIV(Fout[m3],4);

        C_MUL(scratch[0],Fout[m] , *tw1 );
        C_MUL(scratch[1],Fout[m2] , *tw2 );
        C_MUL(scratch[2],Fout[m3] , *tw3 );

        C_SUB( scratch[5] , *Fout, scratch[1] );
        C_ADDTO(*Fout, scratch[1]);
        C_ADD( scratch[3] , scratch[0] , scratch[2] );
        C_SUB( scratch[4] , scratch[0] , scratch[2] );
        C_SUB( Fout[m2], *Fout, scratch[3] );
        tw1 += fstride;
        tw2 += fstride*2;
        tw3 += fstride*3;
        C_ADDTO( *Fout , scratch[3] );

        if(st->inverse) {
            Fout[m].r = scratch[5].r - scratch[4].i;
            Fout[m].i = scratch[5].i + scratch[4].r;
            Fout[m3].r = scratch[5].r + scratch[4].i;
            Fout[m3].i = scratch[5].i - scratch[4].r;
        }else{
            Fout[m].r = scratch[5].r + scratch[4].i;
            Fout[m].i = scratch[5].i - scratch[4].r;
            Fout[m3].r = scratch[5].r - scratch[4].i;
            Fout[m3].i = scratch[5].i + scratch[4].r;
        }
        ++Fout;
    }while(--k);
}

static void kf_bfly3(
         kiss_fft_cpx * Fout,
         const size_t fstride,
         const kiss_fft_cfg st,
         size_t m
         )
{
     size_t k=m;
     const size_t m2 = 2*m;
     kiss_fft_cpx *tw1,*tw2;
     kiss_fft_cpx scratch[5];
     kiss_fft_cpx epi3;
     epi3 = st->twiddles[fstride*m];

     tw1=tw2=st->twiddles;

     do{
         C_FIXDIV(*Fout,3); C_FIXDIV(Fout[m],3); C_FIXDIV(Fout[m2],3);

         C_MUL(scratch[1],Fout[m] , *tw1);
         C_MUL(scratch[2],Fout[m2] , *tw2);

         C_ADD(scratch[3],scratch[1],scratch[2]);
         C_SUB(scratch[0],scratch[1],scratch[2]);
         tw1 += fstride;
         tw2 += fstride*2;

         Fout[m].r = Fout->r - HALF_OF(scratch[3].r);
         Fout[m].i = Fout->i - HALF_OF(scratch[3].i);

         C_MULBYSCALAR( scratch[0] , epi3.i );

         C_ADDTO(*Fout,scratch[3]);

         Fout[m2].r = Fout[m].r + scratch[0].i;
         Fout[m2].i = Fout[m].i - scratch[0].r;

         Fout[m].r -= scratch[0].i;
         Fout[m].i += scratch[0].r;

         ++Fout;
     }while(--k);
}

static void kf_bfly5(
        kiss_fft_cpx * Fout,
        const size_t fstride,
        const kiss_fft_cfg st,
        int m
        )
{
    kiss_fft_cpx *Fout0,*Fout1,*Fout2,*Fout3,*Fout4;
    int u;
    kiss_fft_cpx scratch[13];
    kiss_fft_cpx * twiddles = st->twiddles;
    kiss_fft_cpx *tw;
    kiss_fft_cpx ya,yb;
    ya = twiddles[fstride*m];
    yb = twiddles[fstride*2*m];

    Fout0=Fout;
    Fout1=Fout0+m;
    Fout2=Fout0+2*m;
    Fout3=Fout0+3*m;
    Fout4=Fout0+4*m;

    tw=st->twiddles;
    for ( u=0; u<m; ++u ) {
        C_FIXDIV( *Fout0,5); C_FIXDIV( *Fout1,5); C_FIXDIV( *Fout2,5); C_FIXDIV( *Fout3,5); C_FIXDIV( *Fout4,5);
        scratch[0] = *Fout0;

        C_MUL(scratch[1] ,*Fout1, tw[u*fstride]);
        C_MUL(scratch[2] ,*Fout2, tw[2*u*fstride]);
        C_MUL(scratch[3] ,*Fout3, tw[3*u*fstride]);
        C_MUL(scratch[4] ,*Fout4, tw[4*u*fstride]);

        C_ADD( scratch[7],scratch[1],scratch[4]);
        C_SUB( scratch[10],scratch[1],scratch[4]);
        C_ADD( scratch[8],scratch[2],scratch[3]);
        C_SUB( scratch[9],scratch[2],scratch[3]);

        Fout0->r += scratch[7].r + scratch[8].r;
        Fout0->i += scratch[7].i + scratch[8].i;

        scratch[5].r = scratch[0].r + S_MUL(scratch[7].r,ya.r) + S_MUL(scratch[8].r,yb.r);
        scratch[5].i = scratch[0].i + S_MUL(scratch[7].i,ya.r) + S_MUL(scratch[8].i,yb.r);

        scratch[6].r =  S_MUL(scratch[10].i,ya.i) + S_MUL(scratch[9].i,yb.i);
        scratch[6].i = -S_MUL(scratch[10].r,ya.i) - S_MUL(scratch[9].r,yb.i);

        C_SUB(*Fout1,scratch[5],scratch[6]);
        C_ADD(*Fout4,scratch[5],scratch[6]);

        scratch[11].r = scratch[0].r + S_MUL(scratch[7].r,yb.r) + S_MUL(scratch[8].r,ya.r);
        scratch[11].i = scratch[0].i + S_MUL(scratch[7].i,yb.r) + S_MUL(scratch[8].i,ya.r);
        scratch[12].r = - S_MUL(scratch[10].i,yb.i) + S_MUL(scratch[9].i,ya.i);
        scratch[12].i = S_MUL(scratch[10].r,yb.i) - S_MUL(scratch[9].r,ya.i);

        C_ADD(*Fout2,scratch[11],scratch[12]);
        C_SUB(*Fout3,scratch[11],scratch[12]);

        ++Fout0;++Fout1;++Fout2;++Fout3;++Fout4;
    }
}

/* perform the butterfly for one stage of a mixed radix FFT */
static void kf_bfly_generic(
        kiss_fft_cpx * Fout,
        const size_t fstride,
        const kiss_fft_cfg st,
        int m,
        int p
        )
{
    int u,k,q1,q;
    kiss_fft_cpx * twiddles = st->twiddles;
    kiss_fft_cpx t;
    int Norig = st->nfft;

    kiss_fft_cpx * scratch = (kiss_fft_cpx*)KISS_FFT_TMP_ALLOC(sizeof(kiss_fft_cpx)*p);
    if (scratch == NULL){
        KISS_FFT_ERROR("Memory allocation failed.");
        return;
    }

    for ( u=0; u<m; ++u ) {
        k=u;
        for ( q1=0 ; q1<p ; ++q1 ) {
            scratch[q1] = Fout[ k  ];
            C_FIXDIV(scratch[q1],p);
            k += m;
        }

        k=u;
        for ( q1=0 ; q1<p ; ++q1 ) {
            int twidx=0;
            Fout[ k ] = scratch[0];
            for (q=1;q<p;++q ) {
                twidx += fstride * k;
                if (twidx>=Norig) twidx-=Norig;
                C_MUL(t,scratch[q] , twiddles[twidx] );
                C_ADDTO( Fout[ k ] ,t);
            }
            k += m;
        }
    }
    KISS_FFT_TMP_FREE(scratch);
}

static
void kf_work(
        kiss_fft_cpx * Fout,
        const kiss_fft_cpx * f,
        const size_t fstride,
        int in_stride,
        int * factors,
        const kiss_fft_cfg st
        )
{
    kiss_fft_cpx * Fout_beg=Fout;
    const int p=*factors++; /* the radix  */
    const int m=*factors++; /* stage's fft length/p */
    const kiss_fft_cpx * Fout_end = Fout + p*m;

#ifdef _OPENMP
    // use openmp extensions at the
    // top-level (not recursive)
    if (fstride==1 && p<=5 && m!=1)
    {
        int k;

        // execute the p different work units in different threads
#       pragma omp parallel for
        for (k=0;k<p;++k)
            kf_work( Fout +k*m, f+ fstride*in_stride*k,fstride*p,in_stride,factors,st);
        // all threads have joined by this point

        switch (p) {
            case 2: kf_bfly2(Fout,fstride,st,m); break;
            case 3: kf_bfly3(Fout,fstride,st,m); break;
            case 4: kf_bfly4(Fout,fstride,st,m); break;
            case 5: kf_bfly5(Fout,fstride,st,m); break;
            default: kf_bfly_generic(Fout,fstride,st,m,p); break;
        }
        return;
    }
#endif

    if (m==1) {
        do{
            *Fout = *f;
            f += fstride*in_stride;
        }while(++Fout != Fout_end );
    }else{
        do{
            // recursive call:
            // DFT of size m*p performed by doing
            // p instances of smaller DFTs of size m,
            // each one takes a decimated version of the input
            kf_work( Fout , f, fstride*p, in_stride, factors,st);
            f += fstride*in_stride;
        }while( (Fout += m) != Fout_end );
    }

    Fout=Fout_beg;

    // recombine the p smaller DFTs
    switch (p) {
        case 2: kf_bfly2(Fout,fstride,st,m); break;
        case 3: kf_bfly3(Fout,fstride,st,m); break;
        case 4: kf_bfly4(Fout,fstride,st,m); break;
        case 5: kf_bfly5(Fout,fstride,st,m); break;
        default: kf_bfly_generic(Fout,fstride,st,m,p); break;
    }
}

/*  facbuf is populated by p1,m1,p2,m2, ...
    where
    p[i] * m[i] = m[i-1]
    m0 = n                  */
static
void kf_factor(int n,int * facbuf)
{
    int p=4;
    double floor_sqrt;
    floor_sqrt = floor( sqrt((double)n) );

    /*factor out powers of 4, powers of 2, then any remaining primes */
    do {
        while (n % p) {
            switch (p) {
                case 4: p = 2; break;
                case 2: p = 3; break;
                default: p += 2; break;
            }
            if (p > floor_sqrt)
                p = n;          /* no more factors, skip to end */
        }
        n /= p;
        *facbuf++ = p;
        *facbuf++ = n;
    } while (n > 1);
}

/*
 *
 * User-callable function to allocate all necessary storage space for the fft.
 *
 * The return value is a contiguous block of memory, allocated with malloc.  As such,
 * It can be freed with free(), rather than a kiss_fft-specific function.
 * */
kiss_fft_cfg kiss_fft_alloc(int nfft,int inverse_fft,void * mem,size_t * lenmem )
{
    KISS_FFT_ALIGN_CHECK(mem)

    kiss_fft_cfg st=NULL;
    size_t memneeded = KISS_FFT_ALIGN_SIZE_UP(sizeof(struct kiss_fft_state)
        + sizeof(kiss_fft_cpx)*(nfft-1)); /* twiddle factors*/

    if ( lenmem==NULL ) {
        st = ( kiss_fft_cfg)KISS_FFT_MALLOC( memneeded );
    }else{
        if (mem != NULL && *lenmem >= memneeded)
            st = (kiss_fft_cfg)mem;
        *lenmem = memneeded;
    }
    if (st) {
        int i;
        st->nfft=nfft;
        st->inverse = inverse_fft;

        for (i=0;i<nfft;++i) {
            const double pi=3.141592653589793238462643383279502884197169399375105820974944;
            double phase = -2*pi*i / nfft;
            if (st->inverse)
                phase *= -1;
            kf_cexp(st->twiddles+i, phase );
        }

        kf_factor(nfft,st->factors);
    }
    return st;
}

void kiss_fft_stride(kiss_fft_cfg st,const kiss_fft_cpx *fin,kiss_fft_cpx *fout,int in_stride)
{
    if (fin == fout) {
        //NOTE: this is not really an in-place FFT algorithm.
        //It just performs an out-of-place FFT into a temp buffer
        if (fout == NULL){
            KISS_FFT_ERROR("fout buffer NULL.");
        return;
        }

        kiss_fft_cpx * tmpbuf = (kiss_fft_cpx*)KISS_FFT_TMP_ALLOC( sizeof(kiss_fft_cpx)*st->nfft);
        if (tmpbuf == NULL){
            KISS_FFT_ERROR("Memory allocation error.");
        return;
        }

        kf_work(tmpbuf,fin,1,in_stride, st->factors,st);
        memcpy(fout,tmpbuf,sizeof(kiss_fft_cpx)*st->nfft);
        KISS_FFT_TMP_FREE(tmpbuf);
    }else{
        kf_work( fout, fin, 1,in_stride, st->factors,st );
    }
}

void kiss_fft(kiss_fft_cfg cfg,const kiss_fft_cpx *fin,kiss_fft_cpx *fout)
{
    kiss_fft_stride(cfg,fin,fout,1);
}

void kiss_fft_cleanup(void)
{
    // nothing needed any more
}

int kiss_fft_next_fast_size(int n)
{
    while(1) {
        int m=n;
        while ( (m%2) == 0 ) m/=2;
        while ( (m%3) == 0 ) m/=3;
        while ( (m%5) == 0 ) m/=5;
        if (m<=1)
            break; /* n is completely factorable by twos, threes, and fives */
        n++;
    }
    return n;
}


text:
Valid-License-Identifier: BSD-3-Clause
SPDX-URL: https://spdx.org/licenses/BSD-3-Clause.html
Usage-Guide:
  To use the BSD 3-clause "New" or "Revised" License put the following SPDX
  tag/value pair into a comment according to the placement guidelines in
  the licensing rules documentation:
    SPDX-License-Identifier: BSD-3-Clause
License-Text:

Copyright (c) <year> <owner> . All rights reserved.

Redistribution and use in source and binary forms, with or without modification,
are permitted provided that the following conditions are met:

1. Redistributions of source code must retain the above copyright notice,
this list of conditions and the following disclaimer.

2. 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.

3. Neither the name of the copyright holder nor the names of its 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 HOLDER 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.


text:
Valid-License-Identifier: Unlicense
SPDX-URL: https://spdx.org/licenses/Unlicense.html
Usage-Guide:
  To use the Unlicense put the following SPDX tag/value pair into a
  comment according to the placement guidelines in the licensing rules
  documentation:
    SPDX-License-Identifier: Unlicense
License-Text:

This is free and unencumbered software released into the public domain.

Anyone is free to copy, modify, publish, use, compile, sell, or distribute
this software, either in source code form or as a compiled binary, for any
purpose, commercial or non-commercial, and by any means.

In jurisdictions that recognize copyright laws, the author or authors of this
software dedicate any and all copyright interest in the software to the public
domain. We make this dedication for the benefit of the public at large and
to the detriment of our heirs and successors. We intend this dedication to be
an overt act of relinquishment in perpetuity of all present and future rights
to this software under copyright law.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS
FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS
BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH
THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.

For more information, please refer to <http://unlicense.org/>


text:
/*
 *  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 "_kiss_fft_guts.h"

/*
 Some definitions that allow real or complex filtering
*/
#ifdef REAL_FASTFIR
#define MIN_FFT_LEN 2048
#include "kiss_fftr.h"
typedef kiss_fft_scalar kffsamp_t;
typedef kiss_fftr_cfg kfcfg_t;
#define FFT_ALLOC kiss_fftr_alloc
#define FFTFWD kiss_fftr
#define FFTINV kiss_fftri
#else
#define MIN_FFT_LEN 1024
typedef kiss_fft_cpx kffsamp_t;
typedef kiss_fft_cfg kfcfg_t;
#define FFT_ALLOC kiss_fft_alloc
#define FFTFWD kiss_fft
#define FFTINV kiss_fft
#endif

typedef struct kiss_fastfir_state *kiss_fastfir_cfg;

kiss_fastfir_cfg kiss_fastfir_alloc(const kffsamp_t * imp_resp,size_t n_imp_resp,
        size_t * nfft,void * mem,size_t*lenmem);

/* see do_file_filter for usage */
size_t kiss_fastfir( kiss_fastfir_cfg cfg, kffsamp_t * inbuf, kffsamp_t * outbuf, size_t n, size_t *offset);

static int verbose=0;

struct kiss_fastfir_state{
    size_t nfft;
    size_t ngood;
    kfcfg_t fftcfg;
    kfcfg_t ifftcfg;
    kiss_fft_cpx * fir_freq_resp;
    kiss_fft_cpx * freqbuf;
    size_t n_freq_bins;
    kffsamp_t * tmpbuf;
};

kiss_fastfir_cfg kiss_fastfir_alloc(
        const kffsamp_t * imp_resp,size_t n_imp_resp,
        size_t *pnfft, /* if <= 0, an appropriate size will be chosen */
        void * mem,size_t*lenmem)
{
    kiss_fastfir_cfg st = NULL;
    size_t len_fftcfg,len_ifftcfg;
    size_t memneeded = sizeof(struct kiss_fastfir_state);
    char * ptr;
    size_t i;
    size_t nfft=0;
    float scale;
    int n_freq_bins;
    if (pnfft)
        nfft=*pnfft;

    if (nfft<=0) {
        /* determine fft size as next power of two at least 2x
         the impulse response length*/
        i=n_imp_resp-1;
        nfft=2;
        do{
             nfft<<=1;
        }while (i>>=1);
#ifdef MIN_FFT_LEN
        if ( nfft < MIN_FFT_LEN )
            nfft=MIN_FFT_LEN;
#endif
    }
    if (pnfft)
        *pnfft = nfft;

#ifdef REAL_FASTFIR
    n_freq_bins = nfft/2 + 1;
#else
    n_freq_bins = nfft;
#endif
    /*fftcfg*/
    FFT_ALLOC (nfft, 0, NULL, &len_fftcfg);
    memneeded += len_fftcfg;
    /*ifftcfg*/
    FFT_ALLOC (nfft, 1, NULL, &len_ifftcfg);
    memneeded += len_ifftcfg;
    /* tmpbuf */
    memneeded += sizeof(kffsamp_t) * nfft;
    /* fir_freq_resp */
    memneeded += sizeof(kiss_fft_cpx) * n_freq_bins;
    /* freqbuf */
    memneeded += sizeof(kiss_fft_cpx) * n_freq_bins;

    if (lenmem == NULL) {
        st = (kiss_fastfir_cfg) malloc (memneeded);
    } else {
        if (*lenmem >= memneeded)
            st = (kiss_fastfir_cfg) mem;
        *lenmem = memneeded;
    }
    if (!st)
        return NULL;

    st->nfft = nfft;
    st->ngood = nfft - n_imp_resp + 1;
    st->n_freq_bins = n_freq_bins;
    ptr=(char*)(st+1);

    st->fftcfg = (kfcfg_t)ptr;
    ptr += len_fftcfg;

    st->ifftcfg = (kfcfg_t)ptr;
    ptr += len_ifftcfg;

    st->tmpbuf = (kffsamp_t*)ptr;
    ptr += sizeof(kffsamp_t) * nfft;

    st->freqbuf = (kiss_fft_cpx*)ptr;
    ptr += sizeof(kiss_fft_cpx) * n_freq_bins;

    st->fir_freq_resp = (kiss_fft_cpx*)ptr;
    ptr += sizeof(kiss_fft_cpx) * n_freq_bins;

    FFT_ALLOC (nfft,0,st->fftcfg , &len_fftcfg);
    FFT_ALLOC (nfft,1,st->ifftcfg , &len_ifftcfg);

    memset(st->tmpbuf,0,sizeof(kffsamp_t)*nfft);
    /*zero pad in the middle to left-rotate the impulse response
      This puts the scrap samples at the end of the inverse fft'd buffer */
    st->tmpbuf[0] = imp_resp[ n_imp_resp - 1 ];
    for (i=0;i<n_imp_resp - 1; ++i) {
        st->tmpbuf[ nfft - n_imp_resp + 1 + i ] = imp_resp[ i ];
    }

    FFTFWD(st->fftcfg,st->tmpbuf,st->fir_freq_resp);

    /* TODO: this won't work for fixed point */
    scale = 1.0 / st->nfft;

    for ( i=0; i < st->n_freq_bins; ++i ) {
#ifdef USE_SIMD
        st->fir_freq_resp[i].r *= _mm_set1_ps(scale);
        st->fir_freq_resp[i].i *= _mm_set1_ps(scale);
#else
        st->fir_freq_resp[i].r *= scale;
        st->fir_freq_resp[i].i *= scale;
#endif
    }
    return st;
}

static void fastconv1buf(const kiss_fastfir_cfg st,const kffsamp_t * in,kffsamp_t * out)
{
    size_t i;
    /* multiply the frequency response of the input signal by
     that of the fir filter*/
    FFTFWD( st->fftcfg, in , st->freqbuf );
    for ( i=0; i<st->n_freq_bins; ++i ) {
        kiss_fft_cpx tmpsamp;
        C_MUL(tmpsamp,st->freqbuf[i],st->fir_freq_resp[i]);
        st->freqbuf[i] = tmpsamp;
    }

    /* perform the inverse fft*/
    FFTINV(st->ifftcfg,st->freqbuf,out);
}

/* n : the size of inbuf and outbuf in samples
   return value: the number of samples completely processed
   n-retval samples should be copied to the front of the next input buffer */
static size_t kff_nocopy(
        kiss_fastfir_cfg st,
        const kffsamp_t * inbuf,
        kffsamp_t * outbuf,
        size_t n)
{
    size_t norig=n;
    while (n >= st->nfft ) {
        fastconv1buf(st,inbuf,outbuf);
        inbuf += st->ngood;
        outbuf += st->ngood;
        n -= st->ngood;
    }
    return norig - n;
}

static
size_t kff_flush(kiss_fastfir_cfg st,const kffsamp_t * inbuf,kffsamp_t * outbuf,size_t n)
{
    size_t zpad=0,ntmp;

    ntmp = kff_nocopy(st,inbuf,outbuf,n);
    n -= ntmp;
    inbuf += ntmp;
    outbuf += ntmp;

    zpad = st->nfft - n;
    memset(st->tmpbuf,0,sizeof(kffsamp_t)*st->nfft );
    memcpy(st->tmpbuf,inbuf,sizeof(kffsamp_t)*n );

    fastconv1buf(st,st->tmpbuf,st->tmpbuf);

    memcpy(outbuf,st->tmpbuf,sizeof(kffsamp_t)*( st->ngood - zpad ));
    return ntmp + st->ngood - zpad;
}

size_t kiss_fastfir(
        kiss_fastfir_cfg vst,
        kffsamp_t * inbuf,
        kffsamp_t * outbuf,
        size_t n_new,
        size_t *offset)
{
    size_t ntot = n_new + *offset;
    if (n_new==0) {
        return kff_flush(vst,inbuf,outbuf,ntot);
    }else{
        size_t nwritten = kff_nocopy(vst,inbuf,outbuf,ntot);
        *offset = ntot - nwritten;
        /*save the unused or underused samples at the front of the input buffer */
        memcpy( inbuf , inbuf+nwritten , *offset * sizeof(kffsamp_t) );
        return nwritten;
    }
}

#ifdef FAST_FILT_UTIL
#include <unistd.h>
#include <sys/types.h>
#include <sys/mman.h>
#include <assert.h>

static
void direct_file_filter(
        FILE * fin,
        FILE * fout,
        const kffsamp_t * imp_resp,
        size_t n_imp_resp)
{
    size_t nlag = n_imp_resp - 1;

    const kffsamp_t *tmph;
    kffsamp_t *buf, *circbuf;
    kffsamp_t outval;
    size_t nread;
    size_t nbuf;
    size_t oldestlag = 0;
    size_t k, tap;
#ifndef REAL_FASTFIR
    kffsamp_t tmp;
#endif

    nbuf = 4096;
    buf = (kffsamp_t *) malloc ( sizeof (kffsamp_t) * nbuf);
    circbuf = (kffsamp_t *) malloc (sizeof (kffsamp_t) * nlag);
    if (!circbuf || !buf) {
        perror("circbuf allocation");
        exit(1);
    }

    if ( fread (circbuf, sizeof (kffsamp_t), nlag, fin) !=  nlag ) {
        perror ("insufficient data to overcome transient");
        exit (1);
    }

    do {
        nread = fread (buf, sizeof (kffsamp_t), nbuf, fin);
        if (nread <= 0)
            break;

        for (k = 0; k < nread; ++k) {
            tmph = imp_resp+nlag;
#ifdef REAL_FASTFIR
# ifdef USE_SIMD
            outval = _mm_set1_ps(0);
#else
            outval = 0;
#endif
            for (tap = oldestlag; tap < nlag; ++tap)
                outval += circbuf[tap] * *tmph--;
            for (tap = 0; tap < oldestlag; ++tap)
                outval += circbuf[tap] * *tmph--;
            outval += buf[k] * *tmph;
#else
# ifdef USE_SIMD
            outval.r = outval.i = _mm_set1_ps(0);
#else
            outval.r = outval.i = 0;
#endif
            for (tap = oldestlag; tap < nlag; ++tap){
                C_MUL(tmp,circbuf[tap],*tmph);
                --tmph;
                C_ADDTO(outval,tmp);
            }

            for (tap = 0; tap < oldestlag; ++tap) {
                C_MUL(tmp,circbuf[tap],*tmph);
                --tmph;
                C_ADDTO(outval,tmp);
            }
            C_MUL(tmp,buf[k],*tmph);
            C_ADDTO(outval,tmp);
#endif

            circbuf[oldestlag++] = buf[k];
            buf[k] = outval;

            if (oldestlag == nlag)
                oldestlag = 0;
        }

        if (fwrite (buf, sizeof (buf[0]), nread, fout) != nread) {
            perror ("short write");
            exit (1);
        }
    } while (nread);
    free (buf);
    free (circbuf);
}

static
void do_file_filter(
        FILE * fin,
        FILE * fout,
        const kffsamp_t * imp_resp,
        size_t n_imp_resp,
        size_t nfft )
{
    int fdout;
    size_t n_samps_buf;

    kiss_fastfir_cfg cfg;
    kffsamp_t *inbuf,*outbuf;
    int nread,nwrite;
    size_t idx_inbuf;

    fdout = fileno(fout);

    cfg=kiss_fastfir_alloc(imp_resp,n_imp_resp,&nfft,0,0);

    /* use length to minimize buffer shift*/
    n_samps_buf = 8*4096/sizeof(kffsamp_t);
    n_samps_buf = nfft + 4*(nfft-n_imp_resp+1);

    if (verbose) fprintf(stderr,"bufsize=%d\n",(int)(sizeof(kffsamp_t)*n_samps_buf) );

    /*allocate space and initialize pointers */
    inbuf = (kffsamp_t*)malloc(sizeof(kffsamp_t)*n_samps_buf);
    outbuf = (kffsamp_t*)malloc(sizeof(kffsamp_t)*n_samps_buf);

    idx_inbuf=0;
    do{
        /* start reading at inbuf[idx_inbuf] */
        nread = fread( inbuf + idx_inbuf, sizeof(kffsamp_t), n_samps_buf - idx_inbuf,fin );

        /* If nread==0, then this is a flush.
            The total number of samples in input is idx_inbuf + nread . */
        nwrite = kiss_fastfir(cfg, inbuf, outbuf,nread,&idx_inbuf) * sizeof(kffsamp_t);
        /* kiss_fastfir moved any unused samples to the front of inbuf and updated idx_inbuf */

        if ( write(fdout, outbuf, nwrite) != nwrite ) {
            perror("short write");
            exit(1);
        }
    }while ( nread );
    free(cfg);
    free(inbuf);
    free(outbuf);
}

int main(int argc,char**argv)
{
    kffsamp_t * h;
    int use_direct=0;
    size_t nh,nfft=0;
    FILE *fin=stdin;
    FILE *fout=stdout;
    FILE *filtfile=NULL;
    while (1) {
        int c=getopt(argc,argv,"n:h:i:o:vd");
        if (c==-1) break;
        switch (c) {
            case 'v':
                verbose=1;
                break;
            case 'n':
                nfft=atoi(optarg);
                break;
            case 'i':
                fin = fopen(optarg,"rb");
                if (fin==NULL) {
                    perror(optarg);
                    exit(1);
                }
                break;
            case 'o':
                fout = fopen(optarg,"w+b");
                if (fout==NULL) {
                    perror(optarg);
                    exit(1);
                }
                break;
            case 'h':
                filtfile = fopen(optarg,"rb");
                if (filtfile==NULL) {
                    perror(optarg);
                    exit(1);
                }
                break;
            case 'd':
                use_direct=1;
                break;
            case '?':
                     fprintf(stderr,"usage options:\n"
                            "\t-n nfft: fft size to use\n"
                            "\t-d : use direct FIR filtering, not fast convolution\n"
                            "\t-i filename: input file\n"
                            "\t-o filename: output(filtered) file\n"
                            "\t-n nfft: fft size to use\n"
                            "\t-h filename: impulse response\n");
                     exit (1);
            default:fprintf(stderr,"bad %c\n",c);break;
        }
    }
    if (filtfile==NULL) {
        fprintf(stderr,"You must supply the FIR coeffs via -h\n");
        exit(1);
    }
    fseek(filtfile,0,SEEK_END);
    nh = ftell(filtfile) / sizeof(kffsamp_t);
    if (verbose) fprintf(stderr,"%d samples in FIR filter\n",(int)nh);
    h = (kffsamp_t*)malloc(sizeof(kffsamp_t)*nh);
    fseek(filtfile,0,SEEK_SET);
    if (fread(h,sizeof(kffsamp_t),nh,filtfile) != nh)
        fprintf(stderr,"short read on filter file\n");

    fclose(filtfile);

    if (use_direct)
        direct_file_filter( fin, fout, h,nh);
    else
        do_file_filter( fin, fout, h,nh,nfft);

    if (fout!=stdout) fclose(fout);
    if (fin!=stdin) fclose(fin);

    return 0;
}
#endif


text:
/*
 *  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 KISS_FTR_H
#define KISS_FTR_H

#include "kiss_fft.h"
#ifdef __cplusplus
extern "C" {
#endif

/*

 Real optimized version can save about 45% cpu time vs. complex fft of a real seq.

 */

typedef struct kiss_fftr_state *kiss_fftr_cfg;

kiss_fftr_cfg kiss_fftr_alloc(int nfft,int inverse_fft,void * mem, size_t * lenmem);
/*
 nfft must be even

 If you don't care to allocate space, use mem = lenmem = NULL
*/

void kiss_fftr(kiss_fftr_cfg cfg,const kiss_fft_scalar *timedata,kiss_fft_cpx *freqdata);
/*
 input timedata has nfft scalar points
 output freqdata has nfft/2+1 complex points
*/

void kiss_fftri(kiss_fftr_cfg cfg,const kiss_fft_cpx *freqdata,kiss_fft_scalar *timedata);
/*
 input freqdata has  nfft/2+1 complex points
 output timedata has nfft scalar points
*/

#define kiss_fftr_free KISS_FFT_FREE

#ifdef __cplusplus
}
#endif
#endif


text:
/*
 *  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


text:
/*
 *  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 "kiss_fftndr.h"
#include "_kiss_fft_guts.h"
#define MAX(x,y) ( ( (x)<(y) )?(y):(x) )

struct kiss_fftndr_state
{
    int dimReal;
    int dimOther;
    kiss_fftr_cfg cfg_r;
    kiss_fftnd_cfg cfg_nd;
    void * tmpbuf;
};

static int prod(const int *dims, int ndims)
{
    int x=1;
    while (ndims--)
        x *= *dims++;
    return x;
}

kiss_fftndr_cfg kiss_fftndr_alloc(const int *dims,int ndims,int inverse_fft,void*mem,size_t*lenmem)
{
    KISS_FFT_ALIGN_CHECK(mem)

    kiss_fftndr_cfg st = NULL;
    size_t nr=0 , nd=0,ntmp=0;
    int dimReal = dims[ndims-1];
    int dimOther = prod(dims,ndims-1);
    size_t memneeded;
    char * ptr = NULL;

    (void)kiss_fftr_alloc(dimReal,inverse_fft,NULL,&nr);
    (void)kiss_fftnd_alloc(dims,ndims-1,inverse_fft,NULL,&nd);
    ntmp =
        MAX( 2*dimOther , dimReal+2) * sizeof(kiss_fft_scalar)  // freq buffer for one pass
        + dimOther*(dimReal+2) * sizeof(kiss_fft_scalar);  // large enough to hold entire input in case of in-place

    memneeded = KISS_FFT_ALIGN_SIZE_UP(sizeof( struct kiss_fftndr_state )) + KISS_FFT_ALIGN_SIZE_UP(nr) + KISS_FFT_ALIGN_SIZE_UP(nd) + KISS_FFT_ALIGN_SIZE_UP(ntmp);

    if (lenmem==NULL) {
        ptr = (char*) malloc(memneeded);
    }else{
        if (*lenmem >= memneeded)
            ptr = (char *)mem;
        *lenmem = memneeded;
    }
    if (ptr==NULL)
        return NULL;

    st = (kiss_fftndr_cfg) ptr;
    memset( st , 0 , memneeded);
    ptr += KISS_FFT_ALIGN_SIZE_UP(sizeof(struct kiss_fftndr_state));

    st->dimReal = dimReal;
    st->dimOther = dimOther;
    st->cfg_r = kiss_fftr_alloc( dimReal,inverse_fft,ptr,&nr);
    ptr += KISS_FFT_ALIGN_SIZE_UP(nr);
    st->cfg_nd = kiss_fftnd_alloc(dims,ndims-1,inverse_fft, ptr,&nd);
    ptr += KISS_FFT_ALIGN_SIZE_UP(nd);
    st->tmpbuf = ptr;

    return st;
}

void kiss_fftndr(kiss_fftndr_cfg st,const kiss_fft_scalar *timedata,kiss_fft_cpx *freqdata)
{
    int k1,k2;
    int dimReal = st->dimReal;
    int dimOther = st->dimOther;
    int nrbins = dimReal/2+1;

    kiss_fft_cpx * tmp1 = (kiss_fft_cpx*)st->tmpbuf;
    kiss_fft_cpx * tmp2 = tmp1 + MAX(nrbins,dimOther);

    // timedata is N0 x N1 x ... x Nk real

    // take a real chunk of data, fft it and place the output at correct intervals
    for (k1=0;k1<dimOther;++k1) {
        kiss_fftr( st->cfg_r, timedata + k1*dimReal , tmp1 ); // tmp1 now holds nrbins complex points
        for (k2=0;k2<nrbins;++k2)
           tmp2[ k2*dimOther+k1 ] = tmp1[k2];
    }

    for (k2=0;k2<nrbins;++k2) {
        kiss_fftnd(st->cfg_nd, tmp2+k2*dimOther, tmp1);  // tmp1 now holds dimOther complex points
        for (k1=0;k1<dimOther;++k1)
            freqdata[ k1*(nrbins) + k2] = tmp1[k1];
    }
}

void kiss_fftndri(kiss_fftndr_cfg st,const kiss_fft_cpx *freqdata,kiss_fft_scalar *timedata)
{
    int k1,k2;
    int dimReal = st->dimReal;
    int dimOther = st->dimOther;
    int nrbins = dimReal/2+1;
    kiss_fft_cpx * tmp1 = (kiss_fft_cpx*)st->tmpbuf;
    kiss_fft_cpx * tmp2 = tmp1 + MAX(nrbins,dimOther);

    for (k2=0;k2<nrbins;++k2) {
        for (k1=0;k1<dimOther;++k1)
            tmp1[k1] = freqdata[ k1*(nrbins) + k2 ];
        kiss_fftnd(st->cfg_nd, tmp1, tmp2+k2*dimOther);
    }

    for (k1=0;k1<dimOther;++k1) {
        for (k2=0;k2<nrbins;++k2)
            tmp1[k2] = tmp2[ k2*dimOther+k1 ];
        kiss_fftri( st->cfg_r,tmp1,timedata + k1*dimReal);
    }
}


text:
/*
 *  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 <stdlib.h>
#include <math.h>
#include <stdio.h>
#include <string.h>
#include <unistd.h>

#include "kiss_fft.h"
#include "kiss_fftndr.h"

static
void fft_file(FILE * fin,FILE * fout,int nfft,int isinverse)
{
    kiss_fft_cfg st;
    kiss_fft_cpx * buf;
    kiss_fft_cpx * bufout;

    buf = (kiss_fft_cpx*)malloc(sizeof(kiss_fft_cpx) * nfft );
    bufout = (kiss_fft_cpx*)malloc(sizeof(kiss_fft_cpx) * nfft );
    st = kiss_fft_alloc( nfft ,isinverse ,0,0);

    while ( fread( buf , sizeof(kiss_fft_cpx) * nfft ,1, fin ) > 0 ) {
        kiss_fft( st , buf ,bufout);
        fwrite( bufout , sizeof(kiss_fft_cpx) , nfft , fout );
    }
    free(st);
    free(buf);
    free(bufout);
}

static
void fft_filend(FILE * fin,FILE * fout,int *dims,int ndims,int isinverse)
{
    kiss_fftnd_cfg st;
    kiss_fft_cpx *buf;
    int dimprod=1,i;
    for (i=0;i<ndims;++i)
        dimprod *= dims[i];

    buf = (kiss_fft_cpx *) malloc (sizeof (kiss_fft_cpx) * dimprod);
    st = kiss_fftnd_alloc (dims, ndims, isinverse, 0, 0);

    while (fread (buf, sizeof (kiss_fft_cpx) * dimprod, 1, fin) > 0) {
        kiss_fftnd (st, buf, buf);
        fwrite (buf, sizeof (kiss_fft_cpx), dimprod, fout);
    }
    free (st);
    free (buf);
}

static
void fft_filend_real(FILE * fin,FILE * fout,int *dims,int ndims,int isinverse)
{
    int dimprod=1,i;
    kiss_fftndr_cfg st;
    void *ibuf;
    void *obuf;
    int insize,outsize; // size in bytes

    for (i=0;i<ndims;++i)
        dimprod *= dims[i];
    insize = outsize = dimprod;
    int rdim = dims[ndims-1];

    if (isinverse)
        insize = insize*2*(rdim/2+1)/rdim;
    else
        outsize = outsize*2*(rdim/2+1)/rdim;

    ibuf = malloc(insize*sizeof(kiss_fft_scalar));
    obuf = malloc(outsize*sizeof(kiss_fft_scalar));

    st = kiss_fftndr_alloc(dims, ndims, isinverse, 0, 0);

    while ( fread (ibuf, sizeof(kiss_fft_scalar), insize,  fin) > 0) {
        if (isinverse) {
            kiss_fftndri(st,
                    (kiss_fft_cpx*)ibuf,
                    (kiss_fft_scalar*)obuf);
        }else{
            kiss_fftndr(st,
                    (kiss_fft_scalar*)ibuf,
                    (kiss_fft_cpx*)obuf);
        }
        fwrite (obuf, sizeof(kiss_fft_scalar), outsize,fout);
    }
    free(st);
    free(ibuf);
    free(obuf);
}

static
void fft_file_real(FILE * fin,FILE * fout,int nfft,int isinverse)
{
    kiss_fftr_cfg st;
    kiss_fft_scalar * rbuf;
    kiss_fft_cpx * cbuf;

    rbuf = (kiss_fft_scalar*)malloc(sizeof(kiss_fft_scalar) * nfft );
    cbuf = (kiss_fft_cpx*)malloc(sizeof(kiss_fft_cpx) * (nfft/2+1) );
    st = kiss_fftr_alloc( nfft ,isinverse ,0,0);

    if (isinverse==0) {
        while ( fread( rbuf , sizeof(kiss_fft_scalar) * nfft ,1, fin ) > 0 ) {
            kiss_fftr( st , rbuf ,cbuf);
            fwrite( cbuf , sizeof(kiss_fft_cpx) , (nfft/2 + 1) , fout );
        }
    }else{
        while ( fread( cbuf , sizeof(kiss_fft_cpx) * (nfft/2+1) ,1, fin ) > 0 ) {
            kiss_fftri( st , cbuf ,rbuf);
            fwrite( rbuf , sizeof(kiss_fft_scalar) , nfft , fout );
        }
    }
    free(st);
    free(rbuf);
    free(cbuf);
}

static
int get_dims(char * arg,int * dims)
{
    char *p0;
    int ndims=0;

    do{
        p0 = strchr(arg,',');
        if (p0)
            *p0++ = '\0';
        dims[ndims++] = atoi(arg);
//         fprintf(stderr,"dims[%d] = %d\n",ndims-1,dims[ndims-1]);
        arg = p0;
    }while (p0);
    return ndims;
}

int main(int argc,char ** argv)
{
    int isinverse=0;
    int isreal=0;
    FILE *fin=stdin;
    FILE *fout=stdout;
    int ndims=1;
    int dims[32];
    dims[0] = 1024; /*default fft size*/

    while (1) {
        int c=getopt(argc,argv,"n:iR");
        if (c==-1) break;
        switch (c) {
            case 'n':
                ndims = get_dims(optarg,dims);
                break;
            case 'i':isinverse=1;break;
            case 'R':isreal=1;break;
            case '?':
                     fprintf(stderr,"usage options:\n"
                            "\t-n d1[,d2,d3...]: fft dimension(s)\n"
                            "\t-i : inverse\n"
                            "\t-R : real input samples, not complex\n");
                     exit (1);
            default:fprintf(stderr,"bad %c\n",c);break;
        }
    }

    if ( optind < argc ) {
        if (strcmp("-",argv[optind]) !=0)
            fin = fopen(argv[optind],"rb");
        ++optind;
    }

    if ( optind < argc ) {
        if ( strcmp("-",argv[optind]) !=0 )
            fout = fopen(argv[optind],"wb");
        ++optind;
    }

    if (ndims==1) {
        if (isreal)
            fft_file_real(fin,fout,dims[0],isinverse);
        else
            fft_file(fin,fout,dims[0],isinverse);
    }else{
        if (isreal)
            fft_filend_real(fin,fout,dims,ndims,isinverse);
        else
            fft_filend(fin,fout,dims,ndims,isinverse);
    }

    if (fout!=stdout) fclose(fout);
    if (fin!=stdin) fclose(fin);

    return 0;
}


text:
/*
 *  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 KISS_FFTND_H
#define KISS_FFTND_H

#include "kiss_fft.h"

#ifdef __cplusplus
extern "C" {
#endif

typedef struct kiss_fftnd_state * kiss_fftnd_cfg;

kiss_fftnd_cfg  kiss_fftnd_alloc(const int *dims,int ndims,int inverse_fft,void*mem,size_t*lenmem);
void kiss_fftnd(kiss_fftnd_cfg  cfg,const kiss_fft_cpx *fin,kiss_fft_cpx *fout);

#ifdef __cplusplus
}
#endif
#endif


text:
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_*


text:
/*
 *  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


text:
/*
 *  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 "kiss_fftnd.h"
#include "_kiss_fft_guts.h"

struct kiss_fftnd_state{
    int dimprod; /* dimsum would be mighty tasty right now */
    int ndims;
    int *dims;
    kiss_fft_cfg *states; /* cfg states for each dimension */
    kiss_fft_cpx * tmpbuf; /*buffer capable of hold the entire input */
};

kiss_fftnd_cfg kiss_fftnd_alloc(const int *dims,int ndims,int inverse_fft,void*mem,size_t*lenmem)
{
    KISS_FFT_ALIGN_CHECK(mem)

    kiss_fftnd_cfg st = NULL;
    int i;
    int dimprod=1;
    size_t memneeded = KISS_FFT_ALIGN_SIZE_UP(sizeof(struct kiss_fftnd_state));
    char * ptr = NULL;

    for (i=0;i<ndims;++i) {
        size_t sublen=0;
        kiss_fft_alloc (dims[i], inverse_fft, NULL, &sublen);
        memneeded += sublen;   /* st->states[i] */
        dimprod *= dims[i];
    }
    memneeded += KISS_FFT_ALIGN_SIZE_UP(sizeof(int) * ndims);/*  st->dims */
    memneeded += KISS_FFT_ALIGN_SIZE_UP(sizeof(void*) * ndims);/* st->states  */
    memneeded += KISS_FFT_ALIGN_SIZE_UP(sizeof(kiss_fft_cpx) * dimprod); /* st->tmpbuf */

    if (lenmem == NULL) {/* allocate for the caller*/
        ptr = (char *) malloc (memneeded);
    } else { /* initialize supplied buffer if big enough */
        if (*lenmem >= memneeded)
            ptr = (char *) mem;
        *lenmem = memneeded; /*tell caller how big struct is (or would be) */
    }
    if (!ptr)
        return NULL; /*malloc failed or buffer too small */

    st = (kiss_fftnd_cfg) ptr;
    st->dimprod = dimprod;
    st->ndims = ndims;
    ptr += KISS_FFT_ALIGN_SIZE_UP(sizeof(struct kiss_fftnd_state));

    st->states = (kiss_fft_cfg *)ptr;
    ptr += KISS_FFT_ALIGN_SIZE_UP(sizeof(void*) * ndims);

    st->dims = (int*)ptr;
    ptr += KISS_FFT_ALIGN_SIZE_UP(sizeof(int) * ndims);

    st->tmpbuf = (kiss_fft_cpx*)ptr;
    ptr += KISS_FFT_ALIGN_SIZE_UP(sizeof(kiss_fft_cpx) * dimprod);

    for (i=0;i<ndims;++i) {
        size_t len;
        st->dims[i] = dims[i];
        kiss_fft_alloc (st->dims[i], inverse_fft, NULL, &len);
        st->states[i] = kiss_fft_alloc (st->dims[i], inverse_fft, ptr,&len);
        ptr += len;
    }
    /*
Hi there!

If you're looking at this particular code, it probably means you've got a brain-dead bounds checker
that thinks the above code overwrites the end of the array.

It doesn't.

-- Mark

P.S.
The below code might give you some warm fuzzies and help convince you.
       */
    if ( ptr - (char*)st != (int)memneeded ) {
        fprintf(stderr,
                "################################################################################\n"
                "Internal error! Memory allocation miscalculation\n"
                "################################################################################\n"
               );
    }
    return st;
}

/*
 This works by tackling one dimension at a time.

 In effect,
 Each stage starts out by reshaping the matrix into a DixSi 2d matrix.
 A Di-sized fft is taken of each column, transposing the matrix as it goes.

Here's a 3-d example:
Take a 2x3x4 matrix, laid out in memory as a contiguous buffer
 [ [ [ a b c d ] [ e f g h ] [ i j k l ] ]
   [ [ m n o p ] [ q r s t ] [ u v w x ] ] ]

Stage 0 ( D=2): treat the buffer as a 2x12 matrix
   [ [a b ... k l]
     [m n ... w x] ]

   FFT each column with size 2.
   Transpose the matrix at the same time using kiss_fft_stride.

   [ [ a+m a-m ]
     [ b+n b-n]
     ...
     [ k+w k-w ]
     [ l+x l-x ] ]

   Note fft([x y]) == [x+y x-y]

Stage 1 ( D=3) treats the buffer (the output of stage D=2) as an 3x8 matrix,
   [ [ a+m a-m b+n b-n c+o c-o d+p d-p ]
     [ e+q e-q f+r f-r g+s g-s h+t h-t ]
     [ i+u i-u j+v j-v k+w k-w l+x l-x ] ]

   And perform FFTs (size=3) on each of the columns as above, transposing
   the matrix as it goes.  The output of stage 1 is
       (Legend: ap = [ a+m e+q i+u ]
                am = [ a-m e-q i-u ] )

   [ [ sum(ap) fft(ap)[0] fft(ap)[1] ]
     [ sum(am) fft(am)[0] fft(am)[1] ]
     [ sum(bp) fft(bp)[0] fft(bp)[1] ]
     [ sum(bm) fft(bm)[0] fft(bm)[1] ]
     [ sum(cp) fft(cp)[0] fft(cp)[1] ]
     [ sum(cm) fft(cm)[0] fft(cm)[1] ]
     [ sum(dp) fft(dp)[0] fft(dp)[1] ]
     [ sum(dm) fft(dm)[0] fft(dm)[1] ]  ]

Stage 2 ( D=4) treats this buffer as a 4*6 matrix,
   [ [ sum(ap) fft(ap)[0] fft(ap)[1] sum(am) fft(am)[0] fft(am)[1] ]
     [ sum(bp) fft(bp)[0] fft(bp)[1] sum(bm) fft(bm)[0] fft(bm)[1] ]
     [ sum(cp) fft(cp)[0] fft(cp)[1] sum(cm) fft(cm)[0] fft(cm)[1] ]
     [ sum(dp) fft(dp)[0] fft(dp)[1] sum(dm) fft(dm)[0] fft(dm)[1] ]  ]

   Then FFTs each column, transposing as it goes.

   The resulting matrix is the 3d FFT of the 2x3x4 input matrix.

   Note as a sanity check that the first element of the final
   stage's output (DC term) is
   sum( [ sum(ap) sum(bp) sum(cp) sum(dp) ] )
   , i.e. the summation of all 24 input elements.

*/
void kiss_fftnd(kiss_fftnd_cfg st,const kiss_fft_cpx *fin,kiss_fft_cpx *fout)
{
    int i,k;
    const kiss_fft_cpx * bufin=fin;
    kiss_fft_cpx * bufout;

    /*arrange it so the last bufout == fout*/
    if ( st->ndims & 1 ) {
        bufout = fout;
        if (fin==fout) {
            memcpy( st->tmpbuf, fin, sizeof(kiss_fft_cpx) * st->dimprod );
            bufin = st->tmpbuf;
        }
    }else
        bufout = st->tmpbuf;

    for ( k=0; k < st->ndims; ++k) {
        int curdim = st->dims[k];
        int stride = st->dimprod / curdim;

        for ( i=0 ; i<stride ; ++i )
            kiss_fft_stride( st->states[k], bufin+i , bufout+i*curdim, stride );

        /*toggle back and forth between the two buffers*/
        if (bufout == st->tmpbuf){
            bufout = fout;
            bufin = st->tmpbuf;
        }else{
            bufout = st->tmpbuf;
            bufin = fout;
        }
    }
}


text:
/*
 *  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 "kiss_fftr.h"
#include "_kiss_fft_guts.h"

struct kiss_fftr_state{
    kiss_fft_cfg substate;
    kiss_fft_cpx * tmpbuf;
    kiss_fft_cpx * super_twiddles;
#ifdef USE_SIMD
    void * pad;
#endif
};

kiss_fftr_cfg kiss_fftr_alloc(int nfft,int inverse_fft,void * mem,size_t * lenmem)
{
	KISS_FFT_ALIGN_CHECK(mem)

    int i;
    kiss_fftr_cfg st = NULL;
    size_t subsize = 0, memneeded;

    if (nfft & 1) {
        KISS_FFT_ERROR("Real FFT optimization must be even.");
        return NULL;
    }
    nfft >>= 1;

    kiss_fft_alloc (nfft, inverse_fft, NULL, &subsize);
    memneeded = sizeof(struct kiss_fftr_state) + subsize + sizeof(kiss_fft_cpx) * ( nfft * 3 / 2);

    if (lenmem == NULL) {
        st = (kiss_fftr_cfg) KISS_FFT_MALLOC (memneeded);
    } else {
        if (*lenmem >= memneeded)
            st = (kiss_fftr_cfg) mem;
        *lenmem = memneeded;
    }
    if (!st)
        return NULL;

    st->substate = (kiss_fft_cfg) (st + 1); /*just beyond kiss_fftr_state struct */
    st->tmpbuf = (kiss_fft_cpx *) (((char *) st->substate) + subsize);
    st->super_twiddles = st->tmpbuf + nfft;
    kiss_fft_alloc(nfft, inverse_fft, st->substate, &subsize);

    for (i = 0; i < nfft/2; ++i) {
        double phase =
            -3.14159265358979323846264338327 * ((double) (i+1) / nfft + .5);
        if (inverse_fft)
            phase *= -1;
        kf_cexp (st->super_twiddles+i,phase);
    }
    return st;
}

void kiss_fftr(kiss_fftr_cfg st,const kiss_fft_scalar *timedata,kiss_fft_cpx *freqdata)
{
    /* input buffer timedata is stored row-wise */
    int k,ncfft;
    kiss_fft_cpx fpnk,fpk,f1k,f2k,tw,tdc;

    if ( st->substate->inverse) {
        KISS_FFT_ERROR("kiss fft usage error: improper alloc");
        return;/* The caller did not call the correct function */
    }

    ncfft = st->substate->nfft;

    /*perform the parallel fft of two real signals packed in real,imag*/
    kiss_fft( st->substate , (const kiss_fft_cpx*)timedata, st->tmpbuf );
    /* The real part of the DC element of the frequency spectrum in st->tmpbuf
     * contains the sum of the even-numbered elements of the input time sequence
     * The imag part is the sum of the odd-numbered elements
     *
     * The sum of tdc.r and tdc.i is the sum of the input time sequence.
     *      yielding DC of input time sequence
     * The difference of tdc.r - tdc.i is the sum of the input (dot product) [1,-1,1,-1...
     *      yielding Nyquist bin of input time sequence
     */

    tdc.r = st->tmpbuf[0].r;
    tdc.i = st->tmpbuf[0].i;
    C_FIXDIV(tdc,2);
    CHECK_OVERFLOW_OP(tdc.r ,+, tdc.i);
    CHECK_OVERFLOW_OP(tdc.r ,-, tdc.i);
    freqdata[0].r = tdc.r + tdc.i;
    freqdata[ncfft].r = tdc.r - tdc.i;
#ifdef USE_SIMD
    freqdata[ncfft].i = freqdata[0].i = _mm_set1_ps(0);
#else
    freqdata[ncfft].i = freqdata[0].i = 0;
#endif

    for ( k=1;k <= ncfft/2 ; ++k ) {
        fpk    = st->tmpbuf[k];
        fpnk.r =   st->tmpbuf[ncfft-k].r;
        fpnk.i = - st->tmpbuf[ncfft-k].i;
        C_FIXDIV(fpk,2);
        C_FIXDIV(fpnk,2);

        C_ADD( f1k, fpk , fpnk );
        C_SUB( f2k, fpk , fpnk );
        C_MUL( tw , f2k , st->super_twiddles[k-1]);

        freqdata[k].r = HALF_OF(f1k.r + tw.r);
        freqdata[k].i = HALF_OF(f1k.i + tw.i);
        freqdata[ncfft-k].r = HALF_OF(f1k.r - tw.r);
        freqdata[ncfft-k].i = HALF_OF(tw.i - f1k.i);
    }
}

void kiss_fftri(kiss_fftr_cfg st,const kiss_fft_cpx *freqdata,kiss_fft_scalar *timedata)
{
    /* input buffer timedata is stored row-wise */
    int k, ncfft;

    if (st->substate->inverse == 0) {
        KISS_FFT_ERROR("kiss fft usage error: improper alloc");
        return;/* The caller did not call the correct function */
    }

    ncfft = st->substate->nfft;

    st->tmpbuf[0].r = freqdata[0].r + freqdata[ncfft].r;
    st->tmpbuf[0].i = freqdata[0].r - freqdata[ncfft].r;
    C_FIXDIV(st->tmpbuf[0],2);

    for (k = 1; k <= ncfft / 2; ++k) {
        kiss_fft_cpx fk, fnkc, fek, fok, tmp;
        fk = freqdata[k];
        fnkc.r = freqdata[ncfft - k].r;
        fnkc.i = -freqdata[ncfft - k].i;
        C_FIXDIV( fk , 2 );
        C_FIXDIV( fnkc , 2 );

        C_ADD (fek, fk, fnkc);
        C_SUB (tmp, fk, fnkc);
        C_MUL (fok, tmp, st->super_twiddles[k-1]);
        C_ADD (st->tmpbuf[k],     fek, fok);
        C_SUB (st->tmpbuf[ncfft - k], fek, fok);
#ifdef USE_SIMD
        st->tmpbuf[ncfft - k].i *= _mm_set1_ps(-1.0);
#else
        st->tmpbuf[ncfft - k].i *= -1;
#endif
    }
    kiss_fft (st->substate, st->tmpbuf, (kiss_fft_cpx *) timedata);
}


text:
/*
 *  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 KISS_NDR_H
#define KISS_NDR_H

#include "kiss_fft.h"
#include "kiss_fftr.h"
#include "kiss_fftnd.h"

#ifdef __cplusplus
extern "C" {
#endif

typedef struct kiss_fftndr_state *kiss_fftndr_cfg;

kiss_fftndr_cfg  kiss_fftndr_alloc(const int *dims,int ndims,int inverse_fft,void*mem,size_t*lenmem);
/*
 dims[0] must be even

 If you don't care to allocate space, use mem = lenmem = NULL
*/

void kiss_fftndr(
        kiss_fftndr_cfg cfg,
        const kiss_fft_scalar *timedata,
        kiss_fft_cpx *freqdata);
/*
 input timedata has dims[0] X dims[1] X ... X  dims[ndims-1] scalar points
 output freqdata has dims[0] X dims[1] X ... X  dims[ndims-1]/2+1 complex points
*/

void kiss_fftndri(
        kiss_fftndr_cfg cfg,
        const kiss_fft_cpx *freqdata,
        kiss_fft_scalar *timedata);
/*
 input and output dimensions are the exact opposite of kiss_fftndr
*/

#define kiss_fftndr_free free

#ifdef __cplusplus
}
#endif

#endif


text:
#include <inttypes.h>

class HammingWindow
{
private:
    float *m_coefficients;
    int m_window_size;

public:
    HammingWindow(int window_size);
    ~HammingWindow();
    void applyWindow(float *input);
};


text:
#pragma once

#include <stdint.h>
#include "tools/kiss_fftr.h"

class HammingWindow;

class Processor
{
private:
  HammingWindow *m_hamming_window;
  int m_fft_size;
  int m_window_size;
  kiss_fftr_cfg m_cfg;
  kiss_fft_cpx *m_fft_output;

public:
  float *m_energy;
  float *m_fft_input;

  Processor(int window_size);
  void update(int16_t *samples);
};


text:
#include "Processor.h"
#include "HammingWindow.h"
#include "tools/kiss_fftr.h"

Processor::Processor(int window_size)
{
  m_window_size = window_size;
  // work out the FFT size
  m_fft_size = 1;
  while (m_fft_size < window_size)
  {
    m_fft_size *= 2;
  }
  // create the hamming window to apply to the samples
  m_hamming_window = new HammingWindow(m_fft_size);
  // initialise kiss fftr
  m_cfg = kiss_fftr_alloc(m_fft_size, false, 0, 0);

  m_fft_input = static_cast<float *>(malloc(sizeof(float) * m_fft_size));
  for (int i = 0; i < m_fft_size; i++)
  {
    m_fft_input[i] = 0;
  }
  int energy_size = m_fft_size / 2;
  m_fft_output = static_cast<kiss_fft_cpx *>(malloc(sizeof(kiss_fft_cpx) * energy_size));

  m_energy = static_cast<float *>(malloc(sizeof(float) * window_size / 4));
}

void Processor::update(int16_t *samples)
{
  int offset = (m_fft_size - m_window_size) / 2;
  for (int i = 0; i < m_window_size; i++)
  {
    m_fft_input[offset + i] = samples[i] / 30.0f;
  }
  // apply the hamming window
  m_hamming_window->applyWindow(m_fft_input);
  // do the fft
  kiss_fftr(
      m_cfg,
      m_fft_input,
      reinterpret_cast<kiss_fft_cpx *>(m_fft_output));

  for (int i = 0; i < m_window_size / 4; i++)
  {
    const float real = m_fft_output[i].r;
    const float imag = m_fft_output[i].i;
    m_energy[i] = sqrt((real * real) + (imag * imag));
  }
}


text:
#include <stdlib.h>
#include <math.h>
#include "HammingWindow.h"

HammingWindow::HammingWindow(int window_size)
{
    m_window_size = window_size;
    m_coefficients = static_cast<float *>(malloc(sizeof(float) * m_window_size));
    // create the constants for a hamming window
    const float arg = M_PI * 2.0 / window_size;
    for (int i = 0; i < window_size; i++)
    {
        float float_value = 0.5 - (0.5 * cos(arg * (i + 0.5)));
        // Scale it to fixed point and round it.
        m_coefficients[i] = float_value;
    }
}

HammingWindow::~HammingWindow()
{
    free(m_coefficients);
}

void HammingWindow::applyWindow(float *input)
{
    for (int i = 0; i < m_window_size; i++)
    {
        input[i] = input[i] * m_coefficients[i];
    }
}


text:
#include "I2SMEMSSampler.h"
#include "soc/i2s_reg.h"

I2SMEMSSampler::I2SMEMSSampler(
    i2s_port_t i2s_port,
    i2s_pin_config_t &i2s_pins,
    i2s_config_t i2s_config,
    bool fixSPH0645) : I2SSampler(i2s_port, i2s_config)
{
    m_i2sPins = i2s_pins;
    m_fixSPH0645 = fixSPH0645;
}

void I2SMEMSSampler::configureI2S()
{
    if (m_fixSPH0645)
    {
        // FIXES for SPH0645
        REG_SET_BIT(I2S_TIMING_REG(m_i2sPort), BIT(9));
        REG_SET_BIT(I2S_CONF_REG(m_i2sPort), I2S_RX_MSB_SHIFT);
    }

    i2s_set_pin(m_i2sPort, &m_i2sPins);
}

int I2SMEMSSampler::read(int16_t *samples, int count)
{
    // read from i2s
    int32_t *raw_samples = (int32_t *)malloc(sizeof(int32_t) * count);
    size_t bytes_read = 0;
    i2s_read(m_i2sPort, raw_samples, sizeof(int32_t) * count, &bytes_read, portMAX_DELAY);
    int samples_read = bytes_read / sizeof(int32_t);
    for (int i = 0; i < samples_read; i++)
    {
        samples[i] = (raw_samples[i] & 0xFFFFFFF0) >> 14;
    }
    free(raw_samples);
    return samples_read;
}


text:
#include "ADCSampler.h"

ADCSampler::ADCSampler(adc_unit_t adcUnit, adc1_channel_t adcChannel, const i2s_config_t &i2s_config) : I2SSampler(I2S_NUM_0, i2s_config)
{
    m_adcUnit = adcUnit;
    m_adcChannel = adcChannel;
}

void ADCSampler::configureI2S()
{
    //init ADC pad
    i2s_set_adc_mode(m_adcUnit, m_adcChannel);
    // enable the adc
    i2s_adc_enable(m_i2sPort);
}

int ADCSampler::read(int16_t *samples, int count)
{
    // read from i2s
    size_t bytes_read = 0;
    i2s_read(m_i2sPort, samples, sizeof(int16_t) * count, &bytes_read, portMAX_DELAY);
    int samples_read = bytes_read / sizeof(int16_t);
    for (int i = 0; i < samples_read; i++)
    {
        samples[i] = (2048 - (uint16_t(samples[i]) & 0xfff)) * 15;
    }
    return samples_read;
}


text:
#pragma once

#include <freertos/FreeRTOS.h>
#include <driver/i2s.h>

/**
 * Base Class for both the ADC and I2S sampler
 **/
class I2SSampler
{
protected:
    i2s_port_t m_i2sPort = I2S_NUM_0;
    i2s_config_t m_i2s_config;
    virtual void configureI2S() = 0;
    virtual void processI2SData(void *samples, size_t count){
        // nothing to do for the default case
    };

public:
    I2SSampler(i2s_port_t i2sPort, const i2s_config_t &i2sConfig);
    void start();
    virtual int read(int16_t *samples, int count) = 0;
    void stop();
    int sample_rate()
    {
        return m_i2s_config.sample_rate;
    }
};


text:
#pragma once

#include "I2SSampler.h"

class I2SMEMSSampler : public I2SSampler
{
private:
    i2s_pin_config_t m_i2sPins;
    bool m_fixSPH0645;

protected:
    void configureI2S();

public:
    I2SMEMSSampler(
        i2s_port_t i2s_port,
        i2s_pin_config_t &i2s_pins,
        i2s_config_t i2s_config,
        bool fixSPH0645 = false);
    virtual int read(int16_t *samples, int count);
};


text:
#pragma once

#include "I2SSampler.h"

class ADCSampler : public I2SSampler
{
private:
    adc_unit_t m_adcUnit;
    adc1_channel_t m_adcChannel;

protected:
    void configureI2S();

public:
    ADCSampler(adc_unit_t adc_unit, adc1_channel_t adc_channel, const i2s_config_t &i2s_config);
    virtual int read(int16_t *samples, int count);
};


text:
#include "I2SSampler.h"
#include "driver/i2s.h"

I2SSampler::I2SSampler(i2s_port_t i2sPort, const i2s_config_t &i2s_config) : m_i2sPort(i2sPort), m_i2s_config(i2s_config)
{
}

void I2SSampler::start()
{
    //install and start i2s driver
    i2s_driver_install(m_i2sPort, &m_i2s_config, 0, NULL);
    // set up the I2S configuration from the subclass
    configureI2S();
}

void I2SSampler::stop()
{
    // stop the i2S driver
    i2s_driver_uninstall(m_i2sPort);
}


text:
This directory is intended for PlatformIO Unit Testing and project tests.

Unit Testing is a software testing method by which individual units of
source code, sets of one or more MCU program modules together with associated
control data, usage procedures, and operating procedures, are tested to
determine whether they are fit for use. Unit testing finds problems early
in the development cycle.

More information about PlatformIO Unit Testing:
- https://docs.platformio.org/page/plus/unit-testing.html


text:
This directory is intended for project header files.

A header file is a file containing C declarations and macro definitions
to be shared between several project source files. You request the use of a
header file in your project source file (C, C++, etc) located in `src` folder
by including it, with the C preprocessing directive `#include'.

```src/main.c

#include "header.h"

int main (void)
{
 ...
}
```

Including a header file produces the same results as copying the header file
into each source file that needs it. Such copying would be time-consuming
and error-prone. With a header file, the related declarations appear
in only one place. If they need to be changed, they can be changed in one
place, and programs that include the header file will automatically use the
new version when next recompiled. The header file eliminates the labor of
finding and changing all the copies as well as the risk that a failure to
find one copy will result in inconsistencies within a program.

In C, the usual convention is to give header files names that end with `.h'.
It is most portable to use only letters, digits, dashes, and underscores in
header file names, and at most one dot.

Read more about using header files in official GCC documentation:

* Include Syntax
* Include Operation
* Once-Only Headers
* Computed Includes

https://gcc.gnu.org/onlinedocs/cpp/Header-Files.html


text:
# Changelog

All notable changes to this project will be documented in this file.

The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/),
and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).

repo_name: weavery/clarc

repo_language: OCaml

mime_type: text/plain


text:
(lang dune 2.7)
(context default)


text:
# This is free and unencumbered software released into the public domain.

require 'yaml'
require 'active_support/core_ext/hash'  # `gem install activesupport`

def parse_examples(filepath, type)
  examples = {}
  File.open(filepath) do |file|
    name = nil
    file.readline; file.readline  # skip URL
    file.read.split("\n\n").each do |section|
      case section
        when /  \$ clarc/
          examples[name][:ok] += 1 if section.include?('STOP')
          examples[name][:err] += 1 if section.include?('clarc:')
        when /^([^:]+):\s?(.*)/
          examples[name = $1] = {type: type, ok: 0, err: 0, notes: $2.strip}
        else abort "unknown section: #{section}"
      end
    end
  end
  examples
end

def each_feature(&block)
  features = {}
  features.merge!(parse_examples('test/literals.t', 'literal'))
  features.merge!(parse_examples('test/keywords.t', 'keyword'))
  features.merge!(parse_examples('test/operators.t', 'operator'))
  features.merge!(parse_examples('test/functions.t', 'function'))
  #features.each { |k, v| p [k, v] }; exit
  features.keys.sort.each do |feature|
    block.call(feature, features[feature])
  end
end

task default: %w(README.md)

file "README.md" => %w(test/literals.t test/keywords.t test/operators.t test/functions.t) do |t|
  head = File.read(t.name).split("### Supported Clarity features\n", 2).first
  File.open(t.name, 'w') do |file|
    file.puts head
    file.puts "### Supported Clarity features"
    file.puts
    file.puts ["Feature", "Type", "Status", "Notes"].join(' | ')
    file.puts ["-------", "----", "------", "-----"].join(' | ')
    each_feature do |feature, feature_info|
        file.puts [
          "`#{feature}`",
          feature_info[:type],
          feature_info[:notes] == "Not supported." || feature_info[:notes] == "Not implemented yet." ? "❌" :
            (feature_info[:ok] > 0 ? "✅" : "🚧"),
          feature_info[:notes],
        ].join(' | ').strip
    end
    file.puts
    file.puts "**Legend**: ❌ = not supported. 🚧 = work in progress. ✅ = supported."
  end
end


text: Arto Bendiken <arto@bendiken.net>


text:
(lang dune 2.7)
(name clarc)
(version 0.6.0)
(using menhir 2.0)
(cram enable)


text: # Credits


text:
# Clarc

[![Project license](https://img.shields.io/badge/license-Public%20Domain-blue.svg)](https://unlicense.org)
[![Discord](https://img.shields.io/discord/755852964513579099?label=discord)](https://discord.gg/AvHRCDa)

**Clarc** compiles [Clarity] smart contracts into [Ethereum] virtual machine
(EVM) bytecode.

More specifically, the Clarc compiler, called `clarc`, parses `.clar` files and
compiles them into an equivalent EVM bytecode program that runs on the Ethereum
blockchain.

[![Screencast](https://asciinema.org/a/365265.svg)](https://asciinema.org/a/365265)

*Note: Here be dragons. This is a pre-alpha, work-in-progress
project. Assume nothing works, and you may be pleasantly surprised on
occasion.*

## Installation

### Binary Downloads

The latest release binaries for macOS and Linux are available here:

- [clarc-0.5.0-macos.tar.gz](https://github.com/weavery/clarc/releases/download/0.5.0/clarc-0.5.0-macos.tar.gz)

- [clarc-0.5.0-linux.tar.gz](https://github.com/weavery/clarc/releases/download/0.5.0/clarc-0.5.0-linux.tar.gz)

To install, just download and untar the archive and then copy the resulting
binary to `/usr/local/bin`, as follows:

#### macOS

```bash
wget https://github.com/weavery/clarc/releases/download/0.5.0/clarc-0.5.0-macos.tar.gz

tar xf clarc-0.5.0-macos.tar.gz

sudo install clarc-0.5.0-macos /usr/local/bin/clarc
```

#### Linux

```bash
wget https://github.com/weavery/clarc/releases/download/0.5.0/clarc-0.5.0-linux.tar.gz

tar xf clarc-0.5.0-linux.tar.gz

sudo install clarc-0.5.0-linux /usr/local/bin/clarc
```

### Source Code

If you wish to try out the latest and greatest Clarc, you will need to build
it from source code yourself, which entails setting up an [OCaml] development
environment. Reserve at least half an hour of time and
[see further down](#development) in this document for the particulars.

## Usage

To view Clarc's built-in man page that documents all command-line options, run:

```bash
clarc --help
```

![Manpage](https://github.com/weavery/clarc/blob/master/etc/manpage.jpg)

### Compiling to opcode

To compile the Clarity [`counter.clar`] example contract to programmer-readable
output, known as symbolic opcode, run:

```bash
clarc -t opcode counter.clar
```

The previous writes out the compiled program to standard output, which is
helpful during development and debugging.

### Compiling to bytecode

To compile the Clarity [`counter.clar`] example contract to deployable
bytecode, run:

```bash
clarc -t bytecode counter.clar
```

Alternatively, you can specify an output file name in the usual way, with the
target type inferred from the output file extension (`.bin` for bytecode,
in keeping with the Solidity compiler):

```bash
clarc -o counter.bin counter.clar
```

### Deploying the contract

[MyEtherWallet] is the easiest way to deploy the generated contracts:

![Deploy a Contract](https://github.com/weavery/clarc/blob/master/etc/mew-deploy.jpg)

## Examples

### Supported Contracts

The currently tested contracts, deployed to the public [Ropsten] testnet, are:

| Contract ID | Contract Code | Contract Bytecode | Contract ABI |
| :---------- | :------------ | :---------------- | :----------- |
| [0x8a90b1e93020933295b3bd4ce2317062319351d4] | [`counter.clar`] | [`counter.bin`] | [`counter.json`] |
| [0x2e2487c64b1420111e8d66d751f75f69515c5476] | [`kv-store.clar`] | [`kv-store.bin`] | [`kv-store.json`] |
| [0x9a1b29fc432af1e37af03ed2fee00d742ff7372f] | [`panic.clar`] | [`panic.bin`] | [`panic.json`] |

[MyEtherWallet] is the easiest way to interact with these deployed contracts.

[0x8a90b1e93020933295b3bd4ce2317062319351d4]: https://ropsten.etherscan.io/address/0x8a90b1e93020933295b3bd4ce2317062319351d4
[0x2e2487c64b1420111e8d66d751f75f69515c5476]: https://ropsten.etherscan.io/address/0x2e2487c64b1420111e8d66d751f75f69515c5476
[0x9a1b29fc432af1e37af03ed2fee00d742ff7372f]: https://ropsten.etherscan.io/address/0x9a1b29fc432af1e37af03ed2fee00d742ff7372f

### Example: Counter

#### [`counter.clar`]

```scheme
(define-data-var counter int 0)

(define-read-only (get-counter)
  (ok (var-get counter)))

(define-public (increment)
  (begin
    (var-set counter (+ (var-get counter) 1))
    (ok (var-get counter))))

(define-public (decrement)
  (begin
    (var-set counter (- (var-get counter) 1))
    (ok (var-get counter))))
```

#### `counter.opcode`

```bash
$ clarc counter.clar -t opcode
PUSH1 0x00 PUSH1 0x00 SSTORE PUSH1 0x64 DUP1 PUSH1 0x10 PUSH1 0x00 CODECOPY
PUSH1 0x00 RETURN PUSH1 0xe0 PUSH1 0x02 EXP PUSH1 0x00 CALLDATALOAD DIV DUP1
PUSH4 0x8ada066e EQ PUSH1 0x28 JUMPI DUP1 PUSH4 0xd09de08a EQ PUSH1 0x36
JUMPI DUP1 PUSH4 0x2baeceb7 EQ PUSH1 0x4d JUMPI STOP JUMPDEST POP PUSH1 0x00
SLOAD PUSH1 0x00 MSTORE PUSH1 0x20 PUSH1 0x00 RETURN STOP JUMPDEST POP
PUSH1 0x01 PUSH1 0x00 SLOAD ADD PUSH1 0x00 SSTORE PUSH1 0x00 SLOAD PUSH1 0x00
MSTORE PUSH1 0x20 PUSH1 0x00 RETURN STOP JUMPDEST POP PUSH1 0x01 PUSH1 0x00
SLOAD SUB PUSH1 0x00 SSTORE PUSH1 0x00 SLOAD PUSH1 0x00 MSTORE PUSH1 0x20
PUSH1 0x00 RETURN STOP
```

#### [`counter.bin`]

```bash
$ clarc counter.clar -t bytecode
600060005560648060106000396000f360e060020a6000350480638ada066e146028578063d09de08a1460365780632baeceb714604d57005b5060005460005260206000f3005b5060016000540160005560005460005260206000f3005b5060016000540360005560005460005260206000f300
```

## Design

Clarc is written in [OCaml], an excellent programming language for crafting
compiler toolchains.

Clarc is a standard multi-pass compiler consisting of the following stages:

![Flowchart](https://github.com/weavery/clarc/blob/master/etc/flowchart.png)

The Clarity parser and abstract syntax tree ([AST]) are maintained as a
subproject in an OCaml library called [Clarity.ml]. The library enables anyone
familiar with OCaml to quickly and easily develop more best-of-class tooling
for Clarity contracts.

### Lexical analysis

See Clarity.ml's [`lexer.mll`] for the lexical analyzer source code.

[`lexer.mll`]: https://github.com/weavery/clarity.ml/blob/master/src/lexer.mll

### Syntactic analysis

See Clarity.ml's [`parser.mly`] and [`parse.ml`] for the parser source code.

[`parse.ml`]:   https://github.com/weavery/clarity.ml/blob/master/src/parse.ml
[`parser.mly`]: https://github.com/weavery/clarity.ml/blob/master/src/parser.mly

### Semantic analysis

See Clarity.ml's [`grammar.ml`] for the structure of the Clarity [AST].

[`grammar.ml`]: https://github.com/weavery/clarity.ml/blob/master/src/grammar.ml

## Development

This section documents how to get set up with a development environment for
building Clarc from source code. It is only of interest to people who wish to
contribute to Clarc.

### Prerequisites

The following baseline tooling is required in order to build Clarc from source
code:

- [Git](https://git-scm.com/downloads)

- [OCaml] 4.11+

- [OPAM](https://opam.ocaml.org)

- [Dune](https://dune.build)

- [Docker](https://docs.docker.com/get-docker/) (for release builds only)

We would recommend you *don't* install OCaml from a package manager.

Rather, [get set up with OPAM](https://opam.ocaml.org/doc/Install.html) and
then let OPAM install the correct version of OCaml as follows:

```bash
opam init -c 4.11.1        # if OPAM not yet initialized
opam switch create 4.11.1  # if OPAM already initialized
```

Once OPAM and OCaml are available, install Dune as follows:

```bash
opam install dune
```

### Dependencies

The following OCaml tools and libraries are required in order to build
Clarc from source code:

- [Alcotest](https://opam.ocaml.org/packages/alcotest/)
  for unit tests

- [Clarity.ml] for parsing Clarity code

- [Cmdliner](https://opam.ocaml.org/packages/cmdliner/)
  for the command-line interface

- [Cppo](https://opam.ocaml.org/packages/cppo/)
  for code preprocessing

- [Cryptokit](https://opam.ocaml.org/packages/cryptokit/)
  for the Keccak-256 hash function

- [ISO8601](https://opam.ocaml.org/packages/ISO8601/)
  for date handling

- [Num](https://opam.ocaml.org/packages/num/)
  for 128-bit integers

- [Ocolor](https://opam.ocaml.org/packages/ocolor/)
  for terminal colors

These aforementioned dependencies are all best installed via OPAM:

```bash
opam install -y alcotest cmdliner cppo cryptokit iso8601 num ocolor
opam pin add -y clarity-lang https://github.com/weavery/clarity.ml -k git
```

### Running the program

```bash
alias clarc='dune exec bin/clarc/clarc.exe --'

clarc --help
```

### Installing from source code

```bash
git clone https://github.com/weavery/clarc.git

cd clarc

dune build

sudo install _build/default/bin/clarc/clarc.exe /usr/local/bin/clarc
```

## Acknowledgments

We thank the [Stacks Foundation] for [sponsoring] the development of Clarc.

We thank [Blockstack] and [Algorand] for having developed the Clarity language,
an important evolution for the future of smart contracts.

[Algorand]:          https://algorand.com
[AST]:               https://en.wikipedia.org/wiki/Abstract_syntax_tree
[Blockstack]:        https://blockstack.org
[Clarity]:           https://clarity-lang.org
[Clarity.js]:        https://github.com/weavery/clarity.js
[Clarity.ml]:        https://github.com/weavery/clarity.ml
[Ethereum]:          https://ethereum.org
[IR]:                https://en.wikipedia.org/wiki/Intermediate_representation
[MyEtherWallet]:     https://www.myetherwallet.com/interface/interact-with-contract
[OCaml]:             https://ocaml.org
[Ropsten]:           https://ropsten.etherscan.io
[sponsoring]:        https://github.com/stacksgov/Stacks-Grants/issues/16
[Stacks Foundation]: https://stacks.org

[`counter.clar`]:    https://github.com/weavery/clarc/blob/master/etc/examples/counter.clar
[`counter.bin`]:     https://gist.github.com/artob/1f08c37a55965ff486e6ca99f1ade00d
[`counter.json`]:    https://github.com/weavery/clarc/blob/master/etc/examples/counter.json

[`kv-store.clar`]:   https://github.com/weavery/clarc/blob/master/etc/examples/kv-store.clar
[`kv-store.bin`]:    https://gist.github.com/artob/b0f176f52d6d538d7b195c2fb7f6058a
[`kv-store.json`]:   https://github.com/weavery/clarc/blob/master/etc/examples/kv-store.json

[`panic.clar`]:      https://github.com/weavery/clarc/blob/master/etc/examples/panic.clar
[`panic.bin`]:       https://gist.github.com/artob/945397b444402f6bea7512993608b02c
[`panic.json`]:      https://github.com/weavery/clarc/blob/master/etc/examples/panic.json

## Status

### Supported Clarity features

Feature | Type | Status | Notes
------- | ---- | ------ | -----
`*` | function | ✅ | For two parameters. Without overflow checking.
`+` | function | ✅ | For two parameters. Without overflow checking.
`-` | function | ✅ | For two parameters. Without underflow checking.
`/` | function | ✅ | For two parameters. Without division-by-zero checking.
`<` | function | ✅ |
`<=` | function | ✅ |
`>` | function | ✅ |
`>=` | function | ✅ |
`and` | operator | ✅ | For two parameters.
`append` | function | ✅ |
`as-contract` | operator | 🚧 |
`as-max-len?` | operator | 🚧 |
`asserts!` | function | ✅ |
`at-block` | operator | ❌ | Not implemented yet.
`begin` | operator | 🚧 |
`block-height` | keyword | ✅ |
`buff` | literal | ✅ |
`burn-block-height` | keyword | ✅ |
`concat` | function | ✅ | Only for lists.
`contract-call?` | operator | ❌ | Not implemented yet.
`contract-caller` | keyword | ✅ |
`contract-of` | operator | ❌ | Not implemented yet.
`default-to` | function | ✅ |
`err` | function | ✅ |
`false` | literal | ✅ |
`filter` | function | 🚧 |
`fold` | function | 🚧 |
`ft-get-balance` | function | 🚧 |
`ft-mint?` | function | 🚧 |
`ft-transfer?` | function | 🚧 |
`get` | operator | 🚧 |
`get-block-info?` | operator | ❌ | Not implemented yet.
`hash160` | function | ✅ |
`if` | operator | ✅ |
`impl-trait` | operator | ❌ | Not implemented yet.
`int` | literal | ✅ |
`is-eq` | function | ✅ | For two parameters.
`is-err` | function | ✅ |
`is-in-regtest` | keyword | ✅ |
`is-none` | function | ✅ |
`is-ok` | function | ✅ |
`is-some` | function | ✅ |
`keccak256` | function | ✅ |
`len` | function | ✅ | Only for literals.
`let` | operator | 🚧 |
`list` | function | ✅ |
`map` | function | 🚧 |
`map-delete` | function | 🚧 |
`map-get?` | function | ✅ |
`map-insert` | function | 🚧 |
`map-set` | function | 🚧 |
`match` | operator | ✅ |
`mod` | function | ✅ | Without division-by-zero checking.
`nft-get-owner?` | function | 🚧 |
`nft-mint?` | function | 🚧 |
`nft-transfer?` | function | 🚧 |
`none` | literal | ✅ |
`not` | function | ✅ |
`ok` | function | ✅ |
`or` | operator | ✅ | For two parameters.
`pow` | function | ✅ | Without overflow checking.
`principal` | literal | 🚧 |
`principal-of?` | function | 🚧 |
`print` | function | ✅ | Only for literals. Without a meaningful return value.
`secp256k1-recover?` | function | ❌ | Not implemented yet.
`secp256k1-verify` | function | ❌ | Not implemented yet.
`sha256` | function | ✅ |
`sha512` | function | ❌ | Not implemented yet.
`sha512/256` | function | ❌ | Not implemented yet.
`some` | function | ✅ |
`sqrti` | function | ❌ | Not implemented yet.
`string` | literal | ✅ |
`stx-burn?` | function | ❌ | Not supported.
`stx-get-balance` | function | ❌ | Not supported.
`stx-liquid-supply` | keyword | ❌ | Not supported.
`stx-transfer?` | function | ❌ | Not supported.
`to-int` | function | 🚧 |
`to-uint` | function | 🚧 |
`true` | literal | ✅ |
`try!` | function | ✅ |
`tuple` | operator | 🚧 |
`tx-sender` | keyword | ✅ |
`uint` | literal | ✅ |
`unwrap!` | function | ✅ |
`unwrap-err!` | function | ✅ |
`unwrap-err-panic` | function | ✅ |
`unwrap-panic` | function | ✅ |
`use-trait` | operator | ❌ | Not implemented yet.
`var-get` | operator | ✅ |
`var-set` | operator | ✅ |
`xor` | function | ✅ |

**Legend**: ❌ = not supported. 🚧 = work in progress. ✅ = supported.


text: # To-Dos


text:
# See: https://github.com/weavery/gitstamp-action
---
name: Gitstamp
on:
  push:
    branches:
      - master
jobs:
  gitstamp:
    runs-on: ubuntu-latest
    name: Timestamp commit with Gitstamp
    steps:
      - name: Clone repository
        uses: actions/checkout@v2
      - name: Submit Gitstamp transaction
        uses: weavery/gitstamp-action@v1
        with:
          wallet-key: ${{ secrets.GITSTAMP_KEYFILE }}
          commit-link: true


text:
(executable
 (name clarc)
 (libraries Clar2EVM EVM clarity-lang cmdliner ISO8601 ocolor unix))


text:
(* This is free and unencumbered software released into the public domain. *)

type t =
  | Auto
  | Bytecode
  | Opcode
  | Debug

val of_string : string -> (t, [ `Msg of string ]) result
val to_string : t -> string


text:
(* This is free and unencumbered software released into the public domain. *)

val verbose : bool Cmdliner.Term.t
val files : string list Cmdliner.Term.t
val output : string option Cmdliner.Term.t
val target : Target.t Cmdliner.Term.t
val optimize : int Cmdliner.Term.t
val features : Clar2EVM.Feature.t list Cmdliner.Term.t


text:
(* This is free and unencumbered software released into the public domain. *)

open Cmdliner

let version = "0.6.0"  (* TODO: preprocess from VERSION *)

exception Error of int * string * string

let clarc verbose paths output target _optimize features =
  let _fprintf = Format.fprintf in

  let eprintf = if Unix.isatty (Unix.descr_of_out_channel stderr)
    then Ocolor_format.eprintf
    else Format.eprintf
  in

  let output_channel =
    match output with None -> stdout | Some output_path -> open_out output_path
  in

  let output_formatter = Format.formatter_of_out_channel output_channel in

  let _printf (format : ('a, Format.formatter, unit) format) : 'a =
    if Unix.isatty (Unix.descr_of_out_channel output_channel)
    then
      let formatter' = Ocolor_format.make_formatter output_formatter in
      Format.fprintf (Ocolor_format.unwrap_formatter formatter') format
    else Format.fprintf output_formatter format
  in

  let read_file path =
    let channel = open_in path in
    let contents = really_input_string channel (in_channel_length channel) in
    close_in channel;
    contents
  in

  let guess_target = function
    | None -> None
    | Some output -> begin match Filename.extension output with
        | ".bin" -> Some Target.Bytecode
        | ".opcode" -> Some Target.Opcode
        | _ -> None
      end
  in

  let rec process_program program target =
    match target with
    | Target.Auto -> begin match guess_target output with
        | Some target -> process_program program target
        | None -> process_program program Bytecode
      end
    | Bytecode ->
      let (deployer, program) = Clar2EVM.compile_contract program ~features in
      let program = deployer @ program in
      let printf = Format.fprintf output_formatter in
      printf "@[<h>%a@]@." EVM.print_program_as_bytecode program
    | Opcode ->
      let (deployer, program) = Clar2EVM.compile_contract program ~features in
      let program = deployer @ program in
      let printf = Format.fprintf output_formatter in
      printf "@[<hov>%a@]@." EVM.print_program_as_opcode program
    | Debug ->
      let (_, program) = Clar2EVM.compile_contract program ~features in
      let printf = Format.fprintf output_formatter in
      printf "@[<hov>%a@]@." EVM.print_debug program
  in

  let process_file path =
    if verbose then eprintf "@{<yellow>Compiling %s...@}@." path;
    let input = read_file path in
    let program = Clarity.parse_program input in
    process_program program target
  in

  try `Ok (List.iter process_file paths)
  with Error (code, path, error) -> begin
    eprintf "@{<red>error:@} %s: %s@." path error;
    exit code
  end

let command =
  let doc = "compile Clarity contracts for Ethereum" in
  let man = [
    `S Manpage.s_description;
    `P "$(tname) compiles Clarity contracts into Ethereum virtual machine (EVM) bytecode.";
    `S Manpage.s_bugs; `P "Report any bugs at <https://github.com/weavery/clarc/issues>." ]
  in
  Term.(ret (const clarc $ Options.verbose $ Options.files $ Options.output $ Options.target $ Options.optimize $ Options.features)),
  Term.info "clarc" ~version ~doc ~exits:Term.default_exits ~man

let () = Term.(exit @@ eval command)


text:
(* This is free and unencumbered software released into the public domain. *)

open Cmdliner

let verbose =
  let doc = "Be verbose." in
  Arg.(value & flag & info ["v"; "verbose"] ~doc)

let files =
  Arg.(value & pos_all non_dir_file ["/dev/stdin"] & info [] ~docv:"FILE")

let output =
  let doc = "Specify the output file." in
  Arg.(value & opt (some string) None & info ["o"; "output"] ~docv:"OUTPUT" ~doc)

let target =
  let output_format =
    let parse = Target.of_string in
    let print ppf p = Target.to_string p |> Format.fprintf ppf "%s" in
    Arg.conv ~docv:"TARGET" (parse, print)
  in
  let doc = "Specify the output format: `auto', `bytecode', `opcode'." in
  Arg.(value & opt output_format Target.Auto & info ["t"; "target"] ~docv:"TARGET" ~doc)

let optimize =
  let doc = "Specify the optimization level to use." in
  Arg.(value & opt int 0 & info ["O"; "optimize"] ~docv:"LEVEL" ~doc)

let features =
  let open Clar2EVM in
  let feature_flag =
    let parse = Feature.of_string in
    let print ppf p = Feature.to_string p |> Format.fprintf ppf "%s" in
    Arg.conv ~docv:"FLAG" (parse, print)
  in
  let doc = "Specify optional feature flags: `no-deploy', `only-function='." in
  Arg.(value & opt_all feature_flag [Feature.None] & info ["f"; "feature"] ~docv:"FLAG" ~doc)


text:
(* This is free and unencumbered software released into the public domain. *)

type t =
  | Auto
  | Bytecode
  | Opcode
  | Debug

let of_string = function
  | "auto" -> Ok Auto
  | "bytecode" -> Ok Bytecode
  | "opcode" -> Ok Opcode
  | "debug" -> Ok Debug
  | _ -> Error (`Msg "invalid output format")

let to_string = function
  | Auto -> "auto"
  | Bytecode -> "bytecode"
  | Opcode -> "opcode"
  | Debug -> "debug"


text:
(* This is free and unencumbered software released into the public domain. *)

module Feature : sig
  type t =
    | None
    | NoDeploy
    | OnlyFunction of string

  val of_string : string -> (t, [ `Msg of string ]) result
  val to_string : t -> string
end


text:
(* This is free and unencumbered software released into the public domain. *)

let unsupported_function name type' =
  unsupported (Printf.sprintf "(%s %s)" name (Clarity.type_to_string type'))

let unsupported_function2 name type1 type2 =
  let typename1 = Clarity.type_to_string type1 in
  let typename2 = Clarity.type_to_string type2 in
  unsupported (Printf.sprintf "(%s %s %s)" name typename1 typename2)

let rec type_of_expression = function
  | Clarity.Literal lit -> type_of_literal lit
  | Identifier _ -> Clarity.Unit  (* FIXME: unimplemented "type_of_expression for variable bindings" *)

  | TupleExpression _ -> unimplemented "type_of_expression for tuple expressions"
  | ListExpression [] -> List (0, Unit)
  | ListExpression xs -> List (List.length xs, type_of_expression (List.hd xs))
  | SomeExpression expr -> Optional (type_of_expression expr)
  | Ok expr -> Response (type_of_expression expr, Unit)
  | Err expr -> Response (Unit, type_of_expression expr)

  | Add _ | Sub _ | Mul _ | Div _ | Mod _ | Pow _ | Xor _ -> Clarity.Int
  | Ge _ | Gt _ | Le _ | Lt _ -> Bool

  | And _ | Or _ | Not _ -> Bool
  | DefaultTo (default, _) -> type_of_expression default
  | If (_, then', _) -> type_of_expression then'
  | IsEq _ | IsNone _ | IsSome _ | IsErr _ | IsOk _ -> Bool
  | Len _ -> Uint
  | Let (_, body) -> type_of_expression (match last body with Some x -> x | None -> unreachable ())
  | Match (input, (_, ok_expr), (_, err_expr)) ->
    begin match type_of_expression input with
    | Optional _ | Response _ ->
      begin match type_of_expression ok_expr, type_of_expression err_expr with
      | Response _, Response _ -> Response (Unit, Unit)  (* FIXME *)
      | ok_type, err_type when ok_type = err_type -> Response (Unit, Unit)
      | ok_type, err_type -> unsupported_function2 "match" ok_type err_type
      end
    | t -> unsupported_function "match" t
    end
  | ToInt _ -> Int
  | ToUint _ -> Uint
  | Try _ -> unimplemented "type_of_expression for try!"
  | Unwrap (_, _)
  | UnwrapPanic _
  | UnwrapErr (_, _)
  | UnwrapErrPanic _ -> unimplemented "type_of_expression for unwrap forms"
  | VarGet _ -> unimplemented "type_of_expression for var-get"
  | VarSet _ -> Bool

  | Keyword "block-height" -> Uint
  | Keyword "burn-block-height" -> Uint
  | Keyword "contract-caller" -> Principal
  | Keyword "is-in-regtest" -> Bool
  | Keyword "stx-liquid-supply" -> Uint
  | Keyword "tx-sender" -> Principal
  | Keyword id -> unimplemented (Printf.sprintf "type_of_expression for %s" id)

  | FunctionCall ("append", [list; element]) ->
    begin match type_of_expression list, type_of_expression element with
    | List (n, e1), e2 when e1 = e2 -> List (n + 1, e1)
    | t, e -> unsupported_function2 "append" t e
    end

  | FunctionCall ("concat", [list1; list2]) ->
    begin match type_of_expression list1, type_of_expression list2 with
    | List (n1, e1), List (n2, e2) when e1 = e2 -> List (n1 + n2, e1)
    | t1, t2 -> unsupported_function2 "concat" t1 t2
    end

  | FunctionCall ("get", [_; tuple]) -> type_of_expression tuple
  | FunctionCall ("hash160", _) -> Buff 20
  | FunctionCall ("keccak256", _) -> Buff 32
  | FunctionCall ("map-set", _) -> Bool
  | FunctionCall ("map-get?", _) -> Optional (Tuple [])
  | FunctionCall ("print", [expr]) -> type_of_expression expr
  | FunctionCall ("sha256", _) -> Buff 32
  | FunctionCall ("sha512", _) -> Buff 64
  | FunctionCall ("sha512/256", _) -> Buff 32
  | FunctionCall (id, _) -> unimplemented (Printf.sprintf "type_of_expression for '%s'" id)

and type_of_literal = function
  | Clarity.NoneLiteral -> Clarity.Optional Unit
  | BoolLiteral _ -> Bool
  | IntLiteral _ -> Int
  | UintLiteral _ -> Uint
  | TupleLiteral kvs -> Tuple (List.map (fun (id, lit) -> (id, type_of_literal lit)) kvs)
  | BuffLiteral s -> Buff (String.length s)
  | StringLiteral s -> String (String.length s, Clarity.UTF8)


text:
(library
 (name Clar2EVM)
 (modules Clar2EVM)
 (wrapped false)
 (preprocess (action (run %{bin:cppo} -V OCAML:%{ocaml_version} %{input-file})))
 (preprocessor_deps
  compile.ml
  compile.mli
  features.ml
  features.mli
  typecheck.ml
  typecheck.mli
  utility.ml
  utility.mli)
 (libraries EVM clarity-lang str))


text:
(* This is free and unencumbered software released into the public domain. *)

#include "features.ml"
#include "utility.ml"
#include "typecheck.ml"
#include "compile.ml"


text:
(* This is free and unencumbered software released into the public domain. *)

val type_of_expression : Clarity.expression -> Clarity.type'

val type_of_literal : Clarity.literal -> Clarity.type'


text:
(* This is free and unencumbered software released into the public domain. *)

let unreachable () = failwith "unreachable"

let unimplemented what =
  let message =
    if what = "" then "not implemented yet"
    else Printf.sprintf "%s not implemented yet" what
  in
  failwith message

let unsupported what =
  let message =
    if what = "" then "not supported"
    else Printf.sprintf "%s not supported" what
  in
  failwith message

let rec last = function
  | [] -> None
  | [x] -> Some x
  | _ :: tl -> last tl