text:
/* Modified from  https://github.com/bolderflight/MPU9250/tree/master/src @2019-11-28 */

#ifndef SIMPLE_MPU9250_h
#define SIMPLE_MPU9250_h

#include "Arduino.h"
#include "Wire.h"       /* For I2C */

class SimpleMPU9250
{
public:
    enum GyroRange
    {
      GYRO_RANGE_250DPS,
      GYRO_RANGE_500DPS,
      GYRO_RANGE_1000DPS,
      GYRO_RANGE_2000DPS
    };
    enum AccelRange
    {
      ACCEL_RANGE_2G,
      ACCEL_RANGE_4G,
      ACCEL_RANGE_8G,
      ACCEL_RANGE_16G
    };
    enum DlpfBandwidth
    {
      DLPF_BANDWIDTH_184HZ,
      DLPF_BANDWIDTH_92HZ,
      DLPF_BANDWIDTH_41HZ,
      DLPF_BANDWIDTH_20HZ,
      DLPF_BANDWIDTH_10HZ,
      DLPF_BANDWIDTH_5HZ
    };
    enum LpAccelOdr
    {
      LP_ACCEL_ODR_0_24HZ = 0,
      LP_ACCEL_ODR_0_49HZ = 1,
      LP_ACCEL_ODR_0_98HZ = 2,
      LP_ACCEL_ODR_1_95HZ = 3,
      LP_ACCEL_ODR_3_91HZ = 4,
      LP_ACCEL_ODR_7_81HZ = 5,
      LP_ACCEL_ODR_15_63HZ = 6,
      LP_ACCEL_ODR_31_25HZ = 7,
      LP_ACCEL_ODR_62_50HZ = 8,
      LP_ACCEL_ODR_125HZ = 9,
      LP_ACCEL_ODR_250HZ = 10,
      LP_ACCEL_ODR_500HZ = 11
    };

    SimpleMPU9250(TwoWire &_i2c, const uint8_t _address) {
        this->_i2c    = &_i2c;
        this->_address = _address;
    }

    int8_t begin(void) {
        _i2c->begin();
        _i2c->setClock(400000);          /* 400 kHz _i2c bus */

        /* select clock source to gyro */
        if (writeRegister(PWR_MGMNT_1,CLOCK_SEL_PLL) < 0) {
            return -1;
        }
        // enable _i2c master mode
        if(writeRegister(USER_CTRL,I2C_MST_EN) < 0) {
            return -2;
        }
        // set the _i2c bus speed to 400 kHz
        if(writeRegister(I2C_MST_CTRL,I2C_MST_CLK) < 0) {
            return -3;
        }
        // set AK8963 to Power Down
        writeAK8963Register(AK8963_CNTL1,AK8963_PWR_DOWN);
        // reset the MPU9250
        writeRegister(PWR_MGMNT_1,PWR_RESET);
        // wait for MPU-9250 to come back up
        delay(1);
        // reset the AK8963
        writeAK8963Register(AK8963_CNTL2,AK8963_RESET);
        // select clock source to gyro
        if(writeRegister(PWR_MGMNT_1,CLOCK_SEL_PLL) < 0) {
            return -4;
        }
        // check the WHO AM I byte, expected value is 0x71 (decimal 113) or 0x73 (decimal 115)
        if((whoAmI() != 113)&&(whoAmI() != 115)) {
            return -5;
        }
        // enable accelerometer and gyro
        if(writeRegister(PWR_MGMNT_2,SEN_ENABLE) < 0) {
        return -6;
            }
        // setting accel range to 16G as default
        if(writeRegister(ACCEL_CONFIG,ACCEL_FS_SEL_16G) < 0) {
            return -7;
        }
        _accelScale = G * 16.0f/32767.5f; // setting the accel scale to 16G
        _accelRange = ACCEL_RANGE_16G;
        // setting the gyro range to 2000DPS as default
        if(writeRegister(GYRO_CONFIG,GYRO_FS_SEL_2000DPS) < 0) {
            return -8;
        }
        _gyroScale = 2000.0f/32767.5f * _d2r; // setting the gyro scale to 2000DPS
        _gyroRange = GYRO_RANGE_2000DPS;
        // setting bandwidth to 184Hz as default
        if(writeRegister(ACCEL_CONFIG2,ACCEL_DLPF_184) < 0) {
            return -9;
        }
        if(writeRegister(CONFIG,GYRO_DLPF_184) < 0) { // setting gyro bandwidth to 184Hz
            return -10;
        }
        _bandwidth = DLPF_BANDWIDTH_184HZ;
        // setting the sample rate divider to 0 as default
        if(writeRegister(SMPDIV,0x00) < 0) {
            return -11;
        }
        _srd = 0;
        // enable _i2c master mode
        if(writeRegister(USER_CTRL,I2C_MST_EN) < 0) {
            return -12;
        }
        // set the _i2c bus speed to 400 kHz
        if( writeRegister(I2C_MST_CTRL,I2C_MST_CLK) < 0) {
            return -13;
        }
        // check AK8963 WHO AM I register, expected value is 0x48 (decimal 72)
        if( whoAmIAK8963() != 72 ) {
            return -14;
        }
        /* get the magnetometer calibration */
        // set AK8963 to Power Down
        if(writeAK8963Register(AK8963_CNTL1,AK8963_PWR_DOWN) < 0) {
            return -15;
        }
        delay(100); // long wait between AK8963 mode changes
        // set AK8963 to FUSE ROM access
        if(writeAK8963Register(AK8963_CNTL1,AK8963_FUSE_ROM) < 0) {
            return -16;
        }
        delay(100); // long wait between AK8963 mode changes
        // read the AK8963 ASA registers and compute magnetometer scale factors
        readAK8963Registers(AK8963_ASA,3,_buffer);
        _magScaleX = ((((float)_buffer[0]) - 128.0f)/(256.0f) + 1.0f) * 4912.0f / 32760.0f; // micro Tesla
        _magScaleY = ((((float)_buffer[1]) - 128.0f)/(256.0f) + 1.0f) * 4912.0f / 32760.0f; // micro Tesla
        _magScaleZ = ((((float)_buffer[2]) - 128.0f)/(256.0f) + 1.0f) * 4912.0f / 32760.0f; // micro Tesla
        // set AK8963 to Power Down
        if(writeAK8963Register(AK8963_CNTL1,AK8963_PWR_DOWN) < 0) {
            return -17;
        }
        delay(100); // long wait between AK8963 mode changes
        // set AK8963 to 16 bit resolution, 100 Hz update rate
        if(writeAK8963Register(AK8963_CNTL1,AK8963_CNT_MEAS2) < 0) {
            return -18;
        }
        delay(100); // long wait between AK8963 mode changes
        // select clock source to gyro
        if(writeRegister(PWR_MGMNT_1,CLOCK_SEL_PLL) < 0) {
            return -19;
        }
        // instruct the MPU9250 to get 7 bytes of data from the AK8963 at the sample rate
        readAK8963Registers(AK8963_HXL,7,_buffer);
        // estimate gyro bias
        if (calibrateGyro() < 0) {
            return -20;
        }

        setAccelRange(ACCEL_RANGE_2G);
        setGyroRange(GYRO_RANGE_2000DPS);
        setDlpfBandwidth(DLPF_BANDWIDTH_184HZ);
        setSrd(0);
        // successful init, return 1
        return 1;
    }

    int8_t readSensor(void) {
        // grab the data from the MPU9250
        if (readRegisters(ACCEL_OUT, 21, _buffer) < 0) {
            return -1;
        }

        _axcounts = ((((int16_t)_buffer[0]) << 8) | _buffer[1]);
        _aycounts = ((((int16_t)_buffer[2]) << 8) | _buffer[3]);
        _azcounts = ((((int16_t)_buffer[4]) << 8) | _buffer[5]);
        _gxcounts = ((((int16_t)_buffer[8]) << 8) | _buffer[9]);
        _gycounts = ((((int16_t)_buffer[10]) << 8) | _buffer[11]);
        _gzcounts = ((((int16_t)_buffer[12]) << 8) | _buffer[13]);
        _hxcounts = ((((int16_t)_buffer[15]) << 8) | _buffer[14]);
        _hycounts = ((((int16_t)_buffer[17]) << 8) | _buffer[16]);
        _hzcounts = ((((int16_t)_buffer[19]) << 8) | _buffer[18]);

        _ax = (float)(_axcounts) * _accelScale;
        _ay = (float)(_aycounts) * _accelScale;
        _az = (float)(_azcounts) * _accelScale;
        _gx = ((float)(_gxcounts) * _gyroScale) - _gxb;
        _gy = ((float)(_gycounts) * _gyroScale) - _gyb;
        _gz = ((float)(_gzcounts) * _gyroScale) - _gzb;

        /* Transform the magnetomer to align with the Accelerometer and Gyroscope Axis */
        float _hxAK8963 = ((float)_hxcounts) * _magScaleX;
        float _hyAK8963 = ((float)_hycounts) * _magScaleY;
        float _hzAK8963 = ((float)_hzcounts) * _magScaleZ;
        _hx =  _hyAK8963;
        _hy =  _hxAK8963;
        _hz = -_hzAK8963;

        return 1;
    }

    /* estimates the gyro biases */
    int calibrateGyro() {
        // set the range, bandwidth, and srd
        if (setGyroRange(GYRO_RANGE_250DPS) < 0) {
            return -1;
        }
        if (setDlpfBandwidth(DLPF_BANDWIDTH_20HZ) < 0) {
            return -2;
        }
        if (setSrd(19) < 0) {
            return -3;
        }

        // take samples and find bias
        _gxbD = 0;
        _gybD = 0;
        _gzbD = 0;
        for (size_t i=0; i < _numSamples; i++) {
            readSensor();
            _gxbD += (getGyroX_rads() + _gxb)/((double)_numSamples);
            _gybD += (getGyroY_rads() + _gyb)/((double)_numSamples);
            _gzbD += (getGyroZ_rads() + _gzb)/((double)_numSamples);
            delay(20);
        }
        _gxb = (float)_gxbD;
        _gyb = (float)_gybD;
        _gzb = (float)_gzbD;

        // set the range, bandwidth, and srd back to what they were
        if (setGyroRange(_gyroRange) < 0) {
            return -4;
        }
        if (setDlpfBandwidth(_bandwidth) < 0) {
            return -5;
        }
        if (setSrd(_srd) < 0) {
            return -6;
        }
        return 1;
    }

    /* sets the accelerometer full scale range to values other than default */
    int setAccelRange(AccelRange range) {
        switch(range) {
            case ACCEL_RANGE_2G: {
                // setting the accel range to 2G
                if(writeRegister(ACCEL_CONFIG,ACCEL_FS_SEL_2G) < 0) {
                return -1;
                }
                _accelScale = G * 2.0f/32767.5f; // setting the accel scale to 2G
                break;
            }
            case ACCEL_RANGE_4G: {
                // setting the accel range to 4G
                if(writeRegister(ACCEL_CONFIG,ACCEL_FS_SEL_4G) < 0) {
                    return -1;
                }
                _accelScale = G * 4.0f/32767.5f; // setting the accel scale to 4G
                break;
            }
            case ACCEL_RANGE_8G: {
                // setting the accel range to 8G
                if(writeRegister(ACCEL_CONFIG,ACCEL_FS_SEL_8G) < 0) {
                    return -1;
                }
                _accelScale = G * 8.0f/32767.5f; // setting the accel scale to 8G
                break;
            }
            case ACCEL_RANGE_16G: {
                // setting the accel range to 16G
                if(writeRegister(ACCEL_CONFIG,ACCEL_FS_SEL_16G) < 0) {
                    return -1;
                }
                _accelScale = G * 16.0f/32767.5f; // setting the accel scale to 16G
                break;
            }
        }
        _accelRange = range;
        return 1;
    }

    /* sets the gyro full scale range to values other than default */
    int setGyroRange(GyroRange range) {
        switch(range) {
            case GYRO_RANGE_250DPS: {
                // setting the gyro range to 250DPS
                if(writeRegister(GYRO_CONFIG,GYRO_FS_SEL_250DPS) < 0) {
                    return -1;
                }
                _gyroScale = 250.0f/32767.5f * _d2r; // setting the gyro scale to 250DPS
                break;
            }
            case GYRO_RANGE_500DPS: {
                // setting the gyro range to 500DPS
                if(writeRegister(GYRO_CONFIG,GYRO_FS_SEL_500DPS) < 0) {
                    return -1;
                }
                _gyroScale = 500.0f/32767.5f * _d2r; // setting the gyro scale to 500DPS
                break;
            }
            case GYRO_RANGE_1000DPS: {
                // setting the gyro range to 1000DPS
                if(writeRegister(GYRO_CONFIG,GYRO_FS_SEL_1000DPS) < 0) {
                    return -1;
                }
                _gyroScale = 1000.0f/32767.5f * _d2r; // setting the gyro scale to 1000DPS
                break;
            }
            case GYRO_RANGE_2000DPS: {
                // setting the gyro range to 2000DPS
                if(writeRegister(GYRO_CONFIG,GYRO_FS_SEL_2000DPS) < 0) {
                    return -1;
                }
                _gyroScale = 2000.0f/32767.5f * _d2r; // setting the gyro scale to 2000DPS
                break;
            }
        }
        _gyroRange = range;
        return 1;
    }

    /* sets the DLPF bandwidth to values other than default */
    int setDlpfBandwidth(DlpfBandwidth bandwidth) {
        switch(bandwidth) {
            case DLPF_BANDWIDTH_184HZ: {
                if(writeRegister(ACCEL_CONFIG2,ACCEL_DLPF_184) < 0) { // setting accel bandwidth to 184Hz
                    return -1;
                }
                if(writeRegister(CONFIG,GYRO_DLPF_184) < 0) { // setting gyro bandwidth to 184Hz
                    return -2;
                }
                break;
            }
            case DLPF_BANDWIDTH_92HZ: {
                if(writeRegister(ACCEL_CONFIG2,ACCEL_DLPF_92) < 0) { // setting accel bandwidth to 92Hz
                    return -1;
                }
                if(writeRegister(CONFIG,GYRO_DLPF_92) < 0) { // setting gyro bandwidth to 92Hz
                    return -2;
                }
                break;
            }
            case DLPF_BANDWIDTH_41HZ: {
                if(writeRegister(ACCEL_CONFIG2,ACCEL_DLPF_41) < 0) { // setting accel bandwidth to 41Hz
                    return -1;
                }
                if(writeRegister(CONFIG,GYRO_DLPF_41) < 0) { // setting gyro bandwidth to 41Hz
                    return -2;
                }
                break;
            }
            case DLPF_BANDWIDTH_20HZ: {
                if(writeRegister(ACCEL_CONFIG2,ACCEL_DLPF_20) < 0) { // setting accel bandwidth to 20Hz
                    return -1;
                }
                if(writeRegister(CONFIG,GYRO_DLPF_20) < 0) { // setting gyro bandwidth to 20Hz
                    return -2;
                }
                break;
            }
            case DLPF_BANDWIDTH_10HZ: {
                if(writeRegister(ACCEL_CONFIG2,ACCEL_DLPF_10) < 0) { // setting accel bandwidth to 10Hz
                    return -1;
                }
                if(writeRegister(CONFIG,GYRO_DLPF_10) < 0) { // setting gyro bandwidth to 10Hz
                    return -2;
                }
                break;
            }
            case DLPF_BANDWIDTH_5HZ: {
                if(writeRegister(ACCEL_CONFIG2,ACCEL_DLPF_5) < 0) { // setting accel bandwidth to 5Hz
                    return -1;
                }
                if(writeRegister(CONFIG,GYRO_DLPF_5) < 0) { // setting gyro bandwidth to 5Hz
                    return -2;
                }
                break;
            }
        }
        _bandwidth = bandwidth;
        return 1;
    }

    /* sets the sample rate divider to values other than default */
    int setSrd(uint8_t srd) {
        /* setting the sample rate divider to 19 to facilitate setting up magnetometer */
        if(writeRegister(SMPDIV,19) < 0) { // setting the sample rate divider
            return -1;
        }
        if(srd > 9) {
            // set AK8963 to Power Down
            if(writeAK8963Register(AK8963_CNTL1,AK8963_PWR_DOWN) < 0) {
                return -2;
            }
            delay(100); // long wait between AK8963 mode changes
            // set AK8963 to 16 bit resolution, 8 Hz update rate
            if(writeAK8963Register(AK8963_CNTL1,AK8963_CNT_MEAS1) < 0) {
                return -3;
            }
            delay(100); // long wait between AK8963 mode changes
            // instruct the MPU9250 to get 7 bytes of data from the AK8963 at the sample rate
            readAK8963Registers(AK8963_HXL,7,_buffer);
        } else {
            // set AK8963 to Power Down
            if(writeAK8963Register(AK8963_CNTL1,AK8963_PWR_DOWN) < 0) {
                return -2;
            }
            delay(100); // long wait between AK8963 mode changes
            // set AK8963 to 16 bit resolution, 100 Hz update rate
            if(writeAK8963Register(AK8963_CNTL1,AK8963_CNT_MEAS2) < 0) {
                return -3;
            }
            delay(100); // long wait between AK8963 mode changes
            // instruct the MPU9250 to get 7 bytes of data from the AK8963 at the sample rate
            readAK8963Registers(AK8963_HXL,7,_buffer);
        }
        /* setting the sample rate divider */
        if(writeRegister(SMPDIV,srd) < 0) { // setting the sample rate divider
            return -4;
        }
        _srd = srd;
        return 1;
    }

    /* returns the accelerometer measurement in the x direction, m/s/s */
    float getAccelX_mss() {
        return _ax;
    }

    /* returns the accelerometer measurement in the y direction, m/s/s */
    float getAccelY_mss() {
        return _ay;
    }

    /* returns the accelerometer measurement in the z direction, m/s/s */
    float getAccelZ_mss() {
        return _az;
    }

    /* returns the gyroscope measurement in the x direction, rad/s */
    float getGyroX_rads() {
        return _gx;
    }

    /* returns the gyroscope measurement in the y direction, rad/s */
    float getGyroY_rads() {
        return _gy;
    }

    /* returns the gyroscope measurement in the z direction, rad/s */
    float getGyroZ_rads() {
        return _gz;
    }

    /* returns the magnetometer measurement in the x direction, uT */
    float getMagX_uT() {
        return _hx;
    }

    /* returns the magnetometer measurement in the y direction, uT */
    float getMagY_uT() {
        return _hy;
    }

    /* returns the magnetometer measurement in the z direction, uT */
    float getMagZ_uT() {
        return _hz;
    }

private:
    uint8_t _address;
    TwoWire *_i2c;
    size_t _numBytes; // number of bytes received from I2C

    // track success of interacting with sensor
    int _status;
    // buffer for reading from sensor
    uint8_t _buffer[21];
    // data counts
    int16_t _axcounts,_aycounts,_azcounts;
    int16_t _gxcounts,_gycounts,_gzcounts;
    int16_t _hxcounts,_hycounts,_hzcounts;
    int16_t _tcounts;
    // data buffer
    float _ax, _ay, _az;
    float _gx, _gy, _gz;
    float _hx, _hy, _hz;
    float _t;
    // wake on motion
    uint8_t _womThreshold;
    // scale factors
    float _accelScale;
    float _gyroScale;
    float _magScaleX, _magScaleY, _magScaleZ;
    const float _tempScale = 333.87f;
    const float _tempOffset = 21.0f;
    // configuration
    AccelRange _accelRange;
    GyroRange _gyroRange;
    DlpfBandwidth _bandwidth;
    uint8_t _srd;
    // gyro bias estimation
    size_t _numSamples = 100;
    double _gxbD, _gybD, _gzbD;
    float _gxb, _gyb, _gzb;
    // accel bias and scale factor estimation
    double _axbD, _aybD, _azbD;
    float _axmax, _aymax, _azmax;
    float _axmin, _aymin, _azmin;
    float _axb, _ayb, _azb;
    float _axs = 1.0f;
    float _ays = 1.0f;
    float _azs = 1.0f;
    // magnetometer bias and scale factor estimation
    uint16_t _maxCounts = 1000;
    float _deltaThresh = 0.3f;
    uint8_t _coeff = 8;
    uint16_t _counter;
    float _framedelta, _delta;
    float _hxfilt, _hyfilt, _hzfilt;
    float _hxmax, _hymax, _hzmax;
    float _hxmin, _hymin, _hzmin;
    float _hxb, _hyb, _hzb;
    float _hxs = 1.0f;
    float _hys = 1.0f;
    float _hzs = 1.0f;
    float _avgs;

    /* constants */
    const float G = 9.807f;
    const float _d2r = 3.14159265359f/180.0f;
    /* MPU9250 registers */
    const uint8_t ACCEL_OUT = 0x3B;
    const uint8_t GYRO_OUT = 0x43;
    const uint8_t TEMP_OUT = 0x41;
    const uint8_t EXT_SENS_DATA_00 = 0x49;
    const uint8_t ACCEL_CONFIG = 0x1C;
    const uint8_t ACCEL_FS_SEL_2G = 0x00;
    const uint8_t ACCEL_FS_SEL_4G = 0x08;
    const uint8_t ACCEL_FS_SEL_8G = 0x10;
    const uint8_t ACCEL_FS_SEL_16G = 0x18;
    const uint8_t GYRO_CONFIG = 0x1B;
    const uint8_t GYRO_FS_SEL_250DPS = 0x00;
    const uint8_t GYRO_FS_SEL_500DPS = 0x08;
    const uint8_t GYRO_FS_SEL_1000DPS = 0x10;
    const uint8_t GYRO_FS_SEL_2000DPS = 0x18;
    const uint8_t ACCEL_CONFIG2 = 0x1D;
    const uint8_t ACCEL_DLPF_184 = 0x01;
    const uint8_t ACCEL_DLPF_92 = 0x02;
    const uint8_t ACCEL_DLPF_41 = 0x03;
    const uint8_t ACCEL_DLPF_20 = 0x04;
    const uint8_t ACCEL_DLPF_10 = 0x05;
    const uint8_t ACCEL_DLPF_5 = 0x06;
    const uint8_t CONFIG = 0x1A;
    const uint8_t GYRO_DLPF_184 = 0x01;
    const uint8_t GYRO_DLPF_92 = 0x02;
    const uint8_t GYRO_DLPF_41 = 0x03;
    const uint8_t GYRO_DLPF_20 = 0x04;
    const uint8_t GYRO_DLPF_10 = 0x05;
    const uint8_t GYRO_DLPF_5 = 0x06;
    const uint8_t SMPDIV = 0x19;
    const uint8_t INT_PIN_CFG = 0x37;
    const uint8_t INT_ENABLE = 0x38;
    const uint8_t INT_DISABLE = 0x00;
    const uint8_t INT_PULSE_50US = 0x00;
    const uint8_t INT_WOM_EN = 0x40;
    const uint8_t INT_RAW_RDY_EN = 0x01;
    const uint8_t PWR_MGMNT_1 = 0x6B;
    const uint8_t PWR_CYCLE = 0x20;
    const uint8_t PWR_RESET = 0x80;
    const uint8_t CLOCK_SEL_PLL = 0x01;
    const uint8_t PWR_MGMNT_2 = 0x6C;
    const uint8_t SEN_ENABLE = 0x00;
    const uint8_t DIS_GYRO = 0x07;
    const uint8_t USER_CTRL = 0x6A;
    const uint8_t I2C_MST_EN = 0x20;
    const uint8_t I2C_MST_CLK = 0x0D;
    const uint8_t I2C_MST_CTRL = 0x24;
    const uint8_t I2C_SLV0_ADDR = 0x25;
    const uint8_t I2C_SLV0_REG = 0x26;
    const uint8_t I2C_SLV0_DO = 0x63;
    const uint8_t I2C_SLV0_CTRL = 0x27;
    const uint8_t I2C_SLV0_EN = 0x80;
    const uint8_t I2C_READ_FLAG = 0x80;
    const uint8_t MOT_DETECT_CTRL = 0x69;
    const uint8_t ACCEL_INTEL_EN = 0x80;
    const uint8_t ACCEL_INTEL_MODE = 0x40;
    const uint8_t LP_ACCEL_ODR = 0x1E;
    const uint8_t WOM_THR = 0x1F;
    const uint8_t WHO_AM_I = 0x75;
    const uint8_t FIFO_EN = 0x23;
    const uint8_t FIFO_TEMP = 0x80;
    const uint8_t FIFO_GYRO = 0x70;
    const uint8_t FIFO_ACCEL = 0x08;
    const uint8_t FIFO_MAG = 0x01;
    const uint8_t FIFO_COUNT = 0x72;
    const uint8_t FIFO_READ = 0x74;
    /* AK8963 registers */
    const uint8_t AK8963_I2C_ADDR = 0x0C;
    const uint8_t AK8963_HXL = 0x03;
    const uint8_t AK8963_CNTL1 = 0x0A;
    const uint8_t AK8963_PWR_DOWN = 0x00;
    const uint8_t AK8963_CNT_MEAS1 = 0x12;
    const uint8_t AK8963_CNT_MEAS2 = 0x16;
    const uint8_t AK8963_FUSE_ROM = 0x0F;
    const uint8_t AK8963_CNTL2 = 0x0B;
    const uint8_t AK8963_RESET = 0x01;
    const uint8_t AK8963_ASA = 0x10;
    const uint8_t AK8963_WHO_AM_I = 0x00;

    /* writes a byte to MPU9250 register given a register address and data */
    int writeRegister(uint8_t subAddress, uint8_t data) {
        _i2c->beginTransmission(_address); // open the device
        _i2c->write(subAddress); // write the register address
        _i2c->write(data); // write the data
        _i2c->endTransmission();

        delay(10);

        /* read back the register */
        readRegisters(subAddress, 1, _buffer);
        /* check the read back register against the written register */
        if(_buffer[0] == data) {
            return 1;
        }
        else{
            return -1;
        }
    }

    /* reads registers from MPU9250 given a starting register address, number of bytes, and a pointer to store data */
    int readRegisters(uint8_t subAddress, uint8_t count, uint8_t* dest) {
        _i2c->beginTransmission(_address); // open the device
        _i2c->write(subAddress); // specify the starting register address
        _i2c->endTransmission(false);
        _numBytes = _i2c->requestFrom(_address, count); // specify the number of bytes to receive
        if (_numBytes == count) {
            for(uint8_t i = 0; i < count; i++) {
                dest[i] = _i2c->read();
            }
            return 1;
        } else {
            return -1;
        }
    }

    /* writes a register to the AK8963 given a register address and data */
    int writeAK8963Register(uint8_t subAddress, uint8_t data) {
        // set slave 0 to the AK8963 and set for write
        if (writeRegister(I2C_SLV0_ADDR,AK8963_I2C_ADDR) < 0) {
            return -1;
        }
        // set the register to the desired AK8963 sub address
        if (writeRegister(I2C_SLV0_REG,subAddress) < 0) {
            return -2;
        }
        // store the data for write
        if (writeRegister(I2C_SLV0_DO,data) < 0) {
            return -3;
        }
        // enable _i2c and send 1 byte
        if (writeRegister(I2C_SLV0_CTRL,I2C_SLV0_EN | (uint8_t)1) < 0) {
            return -4;
        }
        // read the register and confirm
        if (readAK8963Registers(subAddress,1,_buffer) < 0) {
            return -5;
        }
        if(_buffer[0] == data) {
            return 1;
        } else{
            return -6;
        }
    }

    /* reads registers from the AK8963 */
    int readAK8963Registers(uint8_t subAddress, uint8_t count, uint8_t* dest) {
        // set slave 0 to the AK8963 and set for read
        if (writeRegister(I2C_SLV0_ADDR,AK8963_I2C_ADDR | I2C_READ_FLAG) < 0) {
            return -1;
        }
        // set the register to the desired AK8963 sub address
        if (writeRegister(I2C_SLV0_REG,subAddress) < 0) {
            return -2;
        }
        // enable _i2c and request the bytes
        if (writeRegister(I2C_SLV0_CTRL,I2C_SLV0_EN | count) < 0) {
            return -3;
        }
        delay(1); // takes some time for these registers to fill
        // read the bytes off the MPU9250 EXT_SENS_DATA registers
        _status = readRegisters(EXT_SENS_DATA_00,count,dest);
        return _status;
    }

    /* gets the MPU9250 WHO_AM_I register value, expected to be 0x71 */
    int whoAmI() {
        // read the WHO AM I register
        if (readRegisters(WHO_AM_I,1,_buffer) < 0) {
            return -1;
        }
        // return the register value
        return _buffer[0];
    }

    /* gets the AK8963 WHO_AM_I register value, expected to be 0x48 */
    int whoAmIAK8963() {
        // read the WHO AM I register
        if (readAK8963Registers(AK8963_WHO_AM_I,1,_buffer) < 0) {
            return -1;
        }
        // return the register value
        return _buffer[0];
    }

};

#endif


text:
import easytello as tello
from pyPS4Controller.controller import Controller

drone = tello.Tello()

class MyController(Controller):

    def __init__(self, **kwargs):
        Controller.__init__(self, **kwargs)

    def on_up_arrow_press(self):
        drone.forward(40)

    def on_down_arrow_press(self):
        drone.back(40)

    def on_right_arrow_press(self):
        drone.right(40)

    def on_left_arrow_press(self):
        drone.left(40)

    def on_circle_press(self):
        drone.land()

    def on_x_press(self):
        drone.takeoff()

    def on_R1_press(self):
        drone.up(40)

    def on_L1_press(self):
        drone.down(40)

    def on_square_press(self):
	drone.cw(90)

    def on_triangle_press(self):
	drone.ccw(90)

controller = MyController(interface="/dev/input/js0", connecting_using_ds4drv=False)
# you can start listening before controller is paired, as long as you pair it within the timeout window
controller.listen(timeout=60)

repo_language: Python

mime_type: text/plain


text:
# Control-DJI-Tello-with-Raspberry-Pi-and-PS4-Controller
## Demo of this project
see demo in my Reddit profile: https://www.reddit.com/r/raspberry_pi/comments/rvpqa9/control_dji_tello_drone_with_ps4_controller_and/?utm_source=share&utm_medium=web2x&context=3
## install easytello library on Raspberry Pi
~~~
pip3 install easytello
~~~
## install pyPS4Controller library on Raspberry Pi
~~~
pip3 install pyps4controller
~~~
you can connect PS4 Controller to Raspberry Pi via USB cable or Bluetooth.
![see this link](https://github.com/ArturSpirin/pyPS4Controller)
## connect to DJI Tello with WiFi


text:
# These are supported funding model platforms

github: # Replace with up to 4 GitHub Sponsors-enabled usernames e.g., [user1, user2]
patreon: camenduru
open_collective: # Replace with a single Open Collective username
ko_fi: camenduru
tidelift: # Replace with a single Tidelift platform-name/package-name e.g., npm/babel
community_bridge: # Replace with a single Community Bridge project-name e.g., cloud-foundry
liberapay: # Replace with a single Liberapay username
issuehunt: # Replace with a single IssueHunt username
otechie: # Replace with a single Otechie username
lfx_crowdfunding: # Replace with a single LFX Crowdfunding project-name e.g., cloud-foundry
custom: # Replace with up to 4 custom sponsorship URLs e.g., ['link1', 'link2']

repo_name: camenduru/comfyui-colab

repo_language: Jupyter Notebook


text:
from subprocess import getoutput

class Terminal:
    def __init__(self):
        pass

    @classmethod
    def INPUT_TYPES(s):
        return {"required": {"image": ("IMAGE",), "text": ("STRING", {"multiline": True})}}
    RETURN_TYPES = ("IMAGE",)
    FUNCTION = "execute"
    def execute(self, image, text):
        out = getoutput(f"{text}")
        print(out)
        return (image, )

    CATEGORY = "Tools"

NODE_CLASS_MAPPINGS = {
    "Terminal": Terminal
}


text:
export const defaultGraph = {
	last_node_id: 9,
	last_link_id: 9,
	nodes: [
		{
			id: 7,
			type: "CLIPTextEncode",
			pos: [413, 389],
			size: { 0: 425.27801513671875, 1: 180.6060791015625 },
			flags: {},
			order: 3,
			mode: 0,
			inputs: [{ name: "clip", type: "CLIP", link: 5 }],
			outputs: [{ name: "CONDITIONING", type: "CONDITIONING", links: [6], slot_index: 0 }],
			properties: {},
			widgets_values: ["bad hands"],
		},
		{
			id: 6,
			type: "CLIPTextEncode",
			pos: [415, 186],
			size: { 0: 422.84503173828125, 1: 164.31304931640625 },
			flags: {},
			order: 2,
			mode: 0,
			inputs: [{ name: "clip", type: "CLIP", link: 3 }],
			outputs: [{ name: "CONDITIONING", type: "CONDITIONING", links: [4], slot_index: 0 }],
			properties: {},
			widgets_values: ["masterpiece best quality girl"],
		},
		{
			id: 5,
			type: "EmptyLatentImage",
			pos: [473, 609],
			size: { 0: 315, 1: 106 },
			flags: {},
			order: 1,
			mode: 0,
			outputs: [{ name: "LATENT", type: "LATENT", links: [2], slot_index: 0 }],
			properties: {},
			widgets_values: [512, 512, 1],
		},
		{
			id: 3,
			type: "KSampler",
			pos: [863, 186],
			size: { 0: 315, 1: 262 },
			flags: {},
			order: 4,
			mode: 0,
			inputs: [
				{ name: "model", type: "MODEL", link: 1 },
				{ name: "positive", type: "CONDITIONING", link: 4 },
				{ name: "negative", type: "CONDITIONING", link: 6 },
				{ name: "latent_image", type: "LATENT", link: 2 },
			],
			outputs: [{ name: "LATENT", type: "LATENT", links: [7], slot_index: 0 }],
			properties: {},
			widgets_values: [8566257, true, 20, 8, "euler", "normal", 1],
		},
		{
			id: 8,
			type: "VAEDecode",
			pos: [1209, 188],
			size: { 0: 210, 1: 46 },
			flags: {},
			order: 5,
			mode: 0,
			inputs: [
				{ name: "samples", type: "LATENT", link: 7 },
				{ name: "vae", type: "VAE", link: 8 },
			],
			outputs: [{ name: "IMAGE", type: "IMAGE", links: [9], slot_index: 0 }],
			properties: {},
		},
		{
			id: 9,
			type: "SaveImage",
			pos: [1451, 189],
			size: { 0: 210, 1: 26 },
			flags: {},
			order: 6,
			mode: 0,
			inputs: [{ name: "images", type: "IMAGE", link: 9 }],
			properties: {},
		},
		{
			id: 4,
			type: "CheckpointLoaderSimple",
			pos: [26, 474],
			size: { 0: 315, 1: 98 },
			flags: {},
			order: 0,
			mode: 0,
			outputs: [
				{ name: "MODEL", type: "MODEL", links: [1], slot_index: 0 },
				{ name: "CLIP", type: "CLIP", links: [3, 5], slot_index: 1 },
				{ name: "VAE", type: "VAE", links: [8], slot_index: 2 },
			],
			properties: {},
			widgets_values: ["v1-5-pruned-emaonly.ckpt"],
		},
	],
	links: [
		[1, 4, 0, 3, 0, "MODEL"],
		[2, 5, 0, 3, 3, "LATENT"],
		[3, 4, 1, 6, 0, "CLIP"],
		[4, 6, 0, 3, 1, "CONDITIONING"],
		[5, 4, 1, 7, 0, "CLIP"],
		[6, 7, 0, 3, 2, "CONDITIONING"],
		[7, 3, 0, 8, 0, "LATENT"],
		[8, 4, 2, 8, 1, "VAE"],
		[9, 8, 0, 9, 0, "IMAGE"],
	],
	groups: [],
	config: {},
	extra: {},
	version: 0.4,
};


text:
; PlatformIO Project Configuration File
;
;   Build options: build flags, source filter
;   Upload options: custom upload port, speed and extra flags
;   Library options: dependencies, extra library storages
;   Advanced options: extra scripting
;
; Please visit documentation for the other options and examples
; https://docs.platformio.org/page/projectconf.html

[env:ttgo]
platform = espressif32
board = esp32dev
framework = arduino
monitor_speed = 115200
monitor_filters = esp32_exception_decoder
lib_deps =
  126@3.4.0 ; FastLED
  FS
  SPI
  SPIFFS
  bodmer/TFT_eSPI @ ^2.3.70
build_flags = -Ofast
  -DUSER_SETUP_LOADED=1
  ; setup for the TTGO T4
  -DTFT_WIDTH=135
  -DTFT_HEIGHT=240
  -DST7789_DRIVER=1
  -DCGRAM_OFFSET=1
  -DTFT_MISO=-1
  -DTFT_MOSI=19
  -DTFT_SCLK=18
  -DTFT_CS=5
  -DTFT_DC=16
  -DTFT_RST=23
  -DTFT_BL=4
  -DTFT_BACKLIGHT_ON=HIGH
  -DLOAD_GLCD=1
  -DLOAD_FONT2=1
  -DLOAD_FONT4=1
  -DLOAD_FONT6=1
  -DLOAD_FONT7=1
  -DLOAD_FONT8=1
  -DLOAD_GFXFF=1
  -DSMOOTH_FONT=1
  -DSPI_FREQUENCY=40000000
  -DSPI_READ_FREQUENCY=6000000

repo_language: C

file_name: README

mime_type: text/x-c


text:
# ESP32 Audio Monitor

This is a port of this [project](https://github.com/atomic14/m5stack-core2-audio-monitor) to work with any ESP32 device with a TFT display.

You can watch a video explainer [here (YouTube)](https://www.youtube.com/watch?v=KaJ4b3HJ3NA) which goes into a bit more detail on the audio capabilities of the device.

[![Demo Video](https://img.youtube.com/vi/KaJ4b3HJ3NA/0.jpg)](https://www.youtube.com/watch?v=KaJ4b3HJ3NA)

You'll need to use PlatformIO to build the project.

To clone the repo use:

```
git clone git@github.com:atomic14/esp32-tft-audio-monitor.git
```

Hopefully, the code should be easy to understand.

# Configuration

All the configuration is located in `config.h` please update this to match the pins you are using for your microphone.

If you want to use the built in ADC instead of a MEMS microphone then comment out line 16 and connect your analog microphone to pin 32. You will need some kind of microphone amplifier board such as the `MAX4466` or `MAX9814` from Adafruit.

To modify the code to work with your particular device you will need to change the settings in `platformio.ini` to match your TFT. The easiest way to do this is to look in the `lib/TFT_eSPI/User_Setups` for your device and copy the required pre-processor defines.

If you are feeling lazy then you can just copy the setup for you TFT over the `User_Setup.h` and remove the defines from the build settings.

If you take the lazy option, then remember to remove `-DUSER_SETUP_LOADED=1`.

```
; setup for the TTGO T4
  -DTFT_WIDTH=135
  -DTFT_HEIGHT=240
  -DST7789_DRIVER=1
  -DCGRAM_OFFSET=1
  -DTFT_MISO=-1
  -DTFT_MOSI=19
  -DTFT_SCLK=18
  -DTFT_CS=5
  -DTFT_DC=16
  -DTFT_RST=23
  -DTFT_BL=4
  -DTFT_BACKLIGHT_ON=HIGH
  -DLOAD_GLCD=1
  -DLOAD_FONT2=1
  -DLOAD_FONT4=1
  -DLOAD_FONT6=1
  -DLOAD_FONT7=1
  -DLOAD_FONT8=1
  -DLOAD_GFXFF=1
  -DSMOOTH_FONT=1
  -DSPI_FREQUENCY=40000000
  -DSPI_READ_FREQUENCY=6000000
```

# Contributing

If you want to add some more visualisations then please open up a pull request and contribute some code.


text:
# These are supported funding model platforms

github: [atomic14]
ko_fi: atomic14


text:
#include "config.h"

// i2s config for using the internal ADC
i2s_config_t i2s_adc_config = {
    .mode = (i2s_mode_t)(I2S_MODE_MASTER | I2S_MODE_RX | I2S_MODE_ADC_BUILT_IN),
    .sample_rate = SAMPLE_RATE,
    .bits_per_sample = I2S_BITS_PER_SAMPLE_16BIT,
    .channel_format = I2S_CHANNEL_FMT_ONLY_LEFT,
    .communication_format = I2S_COMM_FORMAT_I2S_LSB,
    .intr_alloc_flags = ESP_INTR_FLAG_LEVEL1,
    .dma_buf_count = 4,
    .dma_buf_len = 1024,
    .use_apll = false,
    .tx_desc_auto_clear = false,
    .fixed_mclk = 0};

// i2s config for reading from I2S
i2s_config_t i2s_mic_Config = {
    .mode = (i2s_mode_t)(I2S_MODE_MASTER | I2S_MODE_RX),
    .sample_rate = SAMPLE_RATE,
    .bits_per_sample = I2S_BITS_PER_SAMPLE_32BIT,
    .channel_format = I2S_CHANNEL_FMT_ONLY_LEFT,
    .communication_format = I2S_COMM_FORMAT_I2S,
    .intr_alloc_flags = ESP_INTR_FLAG_LEVEL1,
    .dma_buf_count = 4,
    .dma_buf_len = 1024,
    .use_apll = false,
    .tx_desc_auto_clear = false,
    .fixed_mclk = 0};

// i2s microphone pins
i2s_pin_config_t i2s_mic_pins = {
    .bck_io_num = I2S_MIC_SERIAL_CLOCK,
    .ws_io_num = I2S_MIC_LEFT_RIGHT_CLOCK,
    .data_out_num = I2S_PIN_NO_CHANGE,
    .data_in_num = I2S_MIC_SERIAL_DATA};


text:
#pragma once

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

// display toggle button
#define GPIO_BUTTON 0

// sample rate for the system
#define SAMPLE_RATE 16000

// approx 30ms of audio @ 16KHz
#define WINDOW_SIZE 512

// are you using an I2S microphone - comment this out if you want to use an analog mic and ADC input
#define USE_I2S_MIC_INPUT

// I2S Microphone Settings
// Which channel is the I2S microphone on? I2S_CHANNEL_FMT_ONLY_LEFT or I2S_CHANNEL_FMT_ONLY_RIGHT
// Generally they will default to LEFT - but you may need to attach the L/R pin to GND
#define I2S_MIC_CHANNEL I2S_CHANNEL_FMT_ONLY_LEFT
// #define I2S_MIC_CHANNEL I2S_CHANNEL_FMT_ONLY_RIGHT
#define I2S_MIC_SERIAL_CLOCK GPIO_NUM_26
#define I2S_MIC_LEFT_RIGHT_CLOCK GPIO_NUM_22
#define I2S_MIC_SERIAL_DATA GPIO_NUM_21

// Analog Microphone Settings - ADC1_CHANNEL_7 is GPIO35
#define ADC_MIC_CHANNEL ADC1_CHANNEL_7

// i2s config for using the internal ADC
extern i2s_config_t i2s_adc_config;
// i2s config for reading from of I2S
extern i2s_config_t i2s_mic_Config;
// i2s microphone pins
extern i2s_pin_config_t i2s_mic_pins;


text:
#include <TFT_eSPI.h>
#include "Application.h"
#include "I2SSampler.h"
#include "I2SMEMSSampler.h"
#include "ADCSampler.h"
#include "UI/UI.h"
#include "Processor.h"
#include "config.h"

// Task to process samples
void processing_task(void *param)
{
  Application *application = (Application *)param;
  // just sit in a loop processing samples as quickly as possible
  while (true)
  {
    application->process_samples();
  }
}

Application::Application(TFT_eSPI &display)
{
  m_window_size = WINDOW_SIZE;
  m_sample_buffer = (int16_t *)malloc(sizeof(int16_t) * WINDOW_SIZE);
  m_ui = new UI(display, m_window_size);
  m_processor = new Processor(m_window_size);
#ifdef USE_I2S_MIC_INPUT
  m_sampler = new I2SMEMSSampler(I2S_NUM_0, i2s_mic_pins, i2s_mic_Config);
#else
  m_sampler = new ADCSampler(ADC_UNIT_1, ADC_MIC_CHANNEL, i2s_adc_config);
#endif
  pinMode(GPIO_BUTTON, INPUT_PULLUP);
}

void Application::begin()
{
  // set up the processing
  TaskHandle_t processing_task_handle;
  xTaskCreatePinnedToCore(processing_task, "Processing Task", 4096, this, 2, &processing_task_handle, 0);

  // start sampling from i2s device
  m_sampler->start();
}

void Application::process_samples()
{
  // grab the samples
  m_sampler->read(m_sample_buffer, WINDOW_SIZE);
  m_processor->update(m_sample_buffer);
  m_ui->update(m_processor->m_fft_input, m_processor->m_energy);
}

void Application::loop()
{
  if (digitalRead(GPIO_BUTTON) == 0)
  {
    m_ui->toggle_display();
    // delay to allow for the touch to finish
    vTaskDelay(pdMS_TO_TICKS(500));
  }
}


text:
#include <Arduino.h>
#include <TFT_eSPI.h>
#include "Application.h"

Application *application;

void setup()
{
  Serial.begin(115200);
  TFT_eSPI *display = new TFT_eSPI();
  display->begin();
  display->setRotation(1);

  application = new Application(*display);
  application->begin();
}

void loop()
{
  // service the application
  application->loop();
}


text:
#pragma once

#include <driver/i2s.h>

class UI;
class Processor;
class I2SSampler;
class TFT_eSPI;

class Application
{
private:
  int m_window_size;
  int16_t *m_sample_buffer;
  UI *m_ui;
  Processor *m_processor;
  I2SSampler *m_sampler;

  void process_samples();

public:
  Application(TFT_eSPI &display);
  void begin();
  void loop();

  friend void processing_task(void *param);
};


text:
#pragma once

#include "Component.h"

class Palette;
class GraphicEqualiser : public Component
{
private:
  Palette *m_palette;
  int m_num_bins;
  float *bar_chart;
  float *bar_chart_peaks;

public:
  GraphicEqualiser(Palette *palette, int x, int y, int width, int height, int num_bins);
  void update(float *mag);
  void _draw(TFT_eSPI &display);
};


text:
#include "Arduino.h"
#include "Bitmap.h"

Bitmap::Bitmap(int width, int height)
{
  this->width = width;
  this->height = height;
  pixels = (uint16_t *)malloc(width * height * sizeof(uint16_t));
  memset(pixels, 0, width * height * sizeof(uint16_t));
  rows = (uint16_t **)malloc(height * sizeof(uint16_t **));
  for (int i = 0; i < height; i++)
  {
    rows[i] = pixels + width * i;
  }
  for (int y = 0; y < height; y++)
  {
    for (int x = 0; x < width; x++)
    {
      rows[y][x] = 0;
    }
  }
}

void Bitmap::scroll_left()
{
  for (int y = 0; y < height; y++)
  {
    memmove(rows[y], rows[y] + 1, this->width * sizeof(uint16_t) - 1);
  }
}


text:
#pragma once

#pragma once

#include <stdint.h>
#include <algorithm>

#undef min

class Palette
{
protected:
  uint16_t colors[256];

public:
  Palette();
  inline uint16_t get_color(int index)
  {
    return colors[std::max(0, std::min(255, index))];
  }
};


text:
#pragma once

class TFT_eSPI;
class Component
{
public:
  int x, y, width, height;
  bool visible;
  Component(int x, int y, int width, int height) : x(x), y(y), width(width), height(height), visible(true)
  {
  }
  virtual void _draw(TFT_eSPI &display) = 0;
  void draw(TFT_eSPI &display)
  {
    if (visible)
    {
      _draw(display);
    }
  }
};


text:
#pragma once

class Palette;
class Waveform;
class GraphicEqualiser;
class Spectrogram;
class TFT_eSPI;

class UI
{
private:
  Palette *m_palette;
  Waveform *m_waveform;
  GraphicEqualiser *m_graphic_equaliser;
  Spectrogram *m_spectrogram;
  TFT_eSPI &m_display;
  TaskHandle_t m_draw_task_handle;

public:
  UI(TFT_eSPI &display, int window_size);
  void toggle_display();
  void update(float *samples, float *fft);
  void draw();
  friend void drawing_task(void *param);
};


text:
#include <Arduino.h>
#include <TFT_eSPI.h>
#include "Waveform.h"

Waveform::Waveform(TFT_eSPI &display, int x, int y, int width, int height, int num_samples) : Component(x, y, width, height)
{
  m_num_samples = num_samples;
  m_samples = static_cast<float *>(malloc(sizeof(float) * num_samples));
}

void Waveform::update(const float *samples)
{
  memcpy(m_samples, samples, sizeof(float) * m_num_samples);
}

void Waveform::_draw(TFT_eSPI &display)
{
  float x = 0;
  float x_step = (float)width / (float)m_num_samples;
  display.fillRect(x, y, width, height, 0);
  for (int i = 4; i < m_num_samples; i += 4)
  {
    display.drawLine(x, height / 2 + m_samples[i - 4] * 3, x + x_step * 4, height / 2 + m_samples[i] * 3, 0xfff);
    x += x_step * 4;
  }
}


text:
#pragma once
#include <stdint.h>

class Bitmap
{
public:
  uint16_t *pixels;
  uint16_t **rows;
  uint16_t width;
  uint16_t height;
  Bitmap(int width, int height);
  void scroll_left();
};


text:
#include <Arduino.h>
#include <TFT_eSPI.h>
#include "UI.h"
#include "UI/Waveform.h"
#include "UI/Palette.h"
#include "UI/GraphicEqualiser.h"
#include "UI/Spectrogram.h"

// Task to process samples
void drawing_task(void *param)
{
  UI *ui = reinterpret_cast<UI *>(param);
  const TickType_t xMaxBlockTime = pdMS_TO_TICKS(1000);
  while (true)
  {
    // wait to be told to redraw
    uint32_t ulNotificationValue = ulTaskNotifyTake(pdTRUE, xMaxBlockTime);
    if (ulNotificationValue != 0)
    {
      ui->draw();
    }
  }
}

UI::UI(TFT_eSPI &display, int window_size) : m_display(display)
{
  Serial.printf("Display is %d x %d\n", display.width(), display.height());
  m_palette = new Palette();
  m_waveform = new Waveform(display, 0, 0, display.width(), display.height(), window_size);
  m_graphic_equaliser = new GraphicEqualiser(m_palette, 0, 0, display.width(), display.height(), window_size);
  m_spectrogram = new Spectrogram(m_palette, 0, 0, display.width(), display.height());
  // start off with the spectrogram hidden
  m_waveform->visible = true;
  m_spectrogram->visible = false;
  m_graphic_equaliser->visible = false;
  // create a drawing task to update our UI
  xTaskCreatePinnedToCore(drawing_task, "Drawing Task", 4096, this, 1, &m_draw_task_handle, 1);
}

void UI::toggle_display()
{
  bool tmp = m_graphic_equaliser->visible;
  m_graphic_equaliser->visible = m_spectrogram->visible;
  m_spectrogram->visible = m_waveform->visible;
  m_waveform->visible = tmp;
}

void UI::update(float *samples, float *fft)
{
  m_waveform->update(samples);
  m_graphic_equaliser->update(fft);
  m_spectrogram->update(fft);
  xTaskNotify(m_draw_task_handle, 1, eIncrement);
}

unsigned long draw_time = 0;
int draw_count = 0;
void UI::draw()
{
  auto start = millis();
  m_spectrogram->draw(m_display);
  m_graphic_equaliser->draw(m_display);
  m_waveform->draw(m_display);
  auto end = millis();
  draw_time += end - start;
  draw_count++;
  if (draw_count == 20)
  {
    Serial.printf("Drawing time %ld\n", draw_time / 20);
    draw_count = 0;
    draw_time = 0;
  }
}


text:
#pragma once

#include "Component.h"

class Bitmap;
class Palette;

class Spectrogram : public Component
{
private:
  Palette *m_palette;

  Bitmap *bitmap;

public:
  Spectrogram(Palette *palette, int x, int y, int width, int height);
  void update(float *magnitudes);
  void _draw(TFT_eSPI &display);
};


text:
#pragma once

#include "Component.h"

class Waveform : public Component
{
private:
  float *m_samples;
  int m_num_samples;

public:
  Waveform(TFT_eSPI &display, int x, int y, int width, int height, int num_samples);
  void update(const float *samples);
  void _draw(TFT_eSPI &display);
};


text:
#include <Arduino.h>
#include <TFT_eSPI.h>
#include "Spectrogram.h"
#include "Bitmap.h"
#include "Palette.h"

Spectrogram::Spectrogram(Palette *palette, int x, int y, int width, int height) : Component(x, y, width, height)
{
  m_palette = palette;
  this->bitmap = new Bitmap(width, height);
}

void Spectrogram::update(float *mag)
{
  bitmap->scroll_left();
  for (int i = 0; i < bitmap->height; i++)
  {
    bitmap->rows[bitmap->height - i - 1][bitmap->width - 1] = m_palette->get_color(mag[i]);
  }
}

void Spectrogram::_draw(TFT_eSPI &display)
{
  display.setSwapBytes(true);
  display.pushImage(x, y, bitmap->width, bitmap->height, bitmap->pixels);
  display.setSwapBytes(false);
}


text:
#include "Palette.h"
#include <FastLED.h>

uint16_t rgb_to_uint16(float r, float g, float b)
{
  int ri = (int)(r * 31.0f);
  int gi = (int)(g * 63.0f);
  int bi = (int)(b * 31.0f);
  // Serial.printf("%d %d %d\n", ri, gi, bi);
  return (ri << 11) | (gi << 5) | bi;
}

DEFINE_GRADIENT_PALETTE(blackToRed){
    0, 0, 0, 50,      //dark blue
    128, 0, 255, 0,   //green
    164, 255, 220, 0, //bright yellow
    192, 255, 128, 0, //bright orange
    255, 255, 0, 0    //full red
};
CRGBPalette256 blackToRedPal = blackToRed;

Palette::Palette()
{
  for (int i = 0; i < 256; ++i)
  {
    CRGB color = ColorFromPalette(blackToRedPal, i, 255, LINEARBLEND);
    colors[i] = rgb_to_uint16(color.r / 255.0f, color.g / 255.0f, color.b / 255.0f);
  }
}


text:
#include <Arduino.h>
#include <algorithm>
#include <TFT_eSPI.h>
#include "Palette.h"
#include "GraphicEqualiser.h"

#undef min

GraphicEqualiser::GraphicEqualiser(Palette *palette, int x, int y, int width, int height, int num_bins) : Component(x, y, width, height)
{
  m_palette = palette;
  m_num_bins = num_bins;
  bar_chart = static_cast<float *>(malloc(sizeof(float) * num_bins));
  for (int i = 0; i < num_bins; i++)
  {
    bar_chart[i] = 0.0f;
  }
  bar_chart_peaks = static_cast<float *>(malloc(sizeof(float) * num_bins));
  for (int i = 0; i < num_bins; i++)
  {
    bar_chart_peaks[i] = 0.0f;
  }
}

void GraphicEqualiser::update(float *mag)
{
  for (int i = 0; i < m_num_bins; i++)
  {
    float m = mag[i];
    if (m > bar_chart[i])
    {
      bar_chart[i] = m;
    }
    else
    {
      bar_chart[i] = 0.7 * bar_chart[i] + 0.3 * m;
    }
    if (m > bar_chart_peaks[i])
    {
      bar_chart_peaks[i] = m;
    }
    else
    {
      bar_chart_peaks[i] = 0.95 * bar_chart_peaks[i] + 0.05 * m;
    }
  }
}

void GraphicEqualiser::_draw(TFT_eSPI &display)
{
  int x = 0;
  int x_step = int(width / (m_num_bins / 16));
  for (int i = 2; i < m_num_bins / 4; i += 4)
  {
    float ave = 0;
    for (int j = 0; j < 4; j++)
    {
      ave += bar_chart[i + j];
    }
    ave /= 4;
    int bar_value = std::min(height, int(0.5f * ave));
    ave = 0;
    for (int j = 0; j < 4; j++)
    {
      ave += bar_chart_peaks[i + j];
    }
    ave /= 4;
    int peak_value = std::min(height, int(0.5f * ave));
    display.fillRect(x, 0, x_step, height - bar_value - 1, 0);
    display.drawLine(x, height - peak_value - 1, x + x_step - 1, height - peak_value - 1, m_palette->get_color(135 + peak_value));
    display.fillRect(x, height - bar_value - 1, x_step - 1, bar_value, m_palette->get_color(135 + bar_value));
    x += x_step;
  }
  display.fillRect(x, 0, width - x, height, 0);
}


text:
{
    // See http://go.microsoft.com/fwlink/?LinkId=827846
    // for the documentation about the extensions.json format
    "recommendations": [
        "platformio.platformio-ide"
    ]
}


text:
This directory is intended for project specific (private) libraries.
PlatformIO will compile them to static libraries and link into executable file.

The source code of each library should be placed in a an own separate directory
("lib/your_library_name/[here are source files]").

For example, see a structure of the following two libraries `Foo` and `Bar`:

|--lib
|  |
|  |--Bar
|  |  |--docs
|  |  |--examples
|  |  |--src
|  |     |- Bar.c
|  |     |- Bar.h
|  |  |- library.json (optional, custom build options, etc) https://docs.platformio.org/page/librarymanager/config.html
|  |
|  |--Foo
|  |  |- Foo.c
|  |  |- Foo.h
|  |
|  |- README --> THIS FILE
|
|- platformio.ini
|--src
   |- main.c

and a contents of `src/main.c`:
```
#include <Foo.h>
#include <Bar.h>

int main (void)
{
  ...
}

```

PlatformIO Library Dependency Finder will find automatically dependent
libraries scanning project source files.

More information about PlatformIO Library Dependency Finder
- https://docs.platformio.org/page/librarymanager/ldf.html


text:
#ifndef KISSFFT_I32_CLASS_HH
#define KISSFFT_I32_CLASS_HH

#include <complex>
#include <utility>
#include <vector>

// TODO1: substitute complex<type> (behaviour not defined for nonfloats), should be faster
// TODO2: use std:: namespace
// TODO3: make unittests for all ffts (c, cpp, i32)

template <typename DType>
struct complex_s
{
    DType real;
    DType imag;
};

class kissfft_i32
{
private:

    using scalar_type = int32_t;
    using cpx_type    = complex<int32_t>;

    scalar_type _scale_factor;
    std::size_t _nfft;
    bool _inverse;
    std::vector<cpx_type> _twiddles;
    std::vector<std::size_t> _stageRadix;
    std::vector<std::size_t> _stageRemainder;

public:

    // scale_factor: upscale twiddle-factors otherwise they lie between 0..1 (out of range for integer) --> fixed point math
    kissfft_i32(const std::size_t nfft, const bool inverse, const double scale_factor = 1024.0)
            : _scale_factor(scalar_type(scale_factor)), _nfft(nfft), _inverse(inverse)
    {
        // fill twiddle factors
        _twiddles.resize(_nfft);
        const double phinc = (_inverse ? 2 : -2) * acos(-1.0) / _nfft;
        for (std::size_t i = 0; i < _nfft; ++i)
        {
            _twiddles[i] = scale_factor * exp(complex<double>(0, i * phinc));
        }
        //factorize
        //start factoring out 4's, then 2's, then 3,5,7,9,...
        std::size_t n = _nfft;
        std::size_t p = 4;
        do
        {
            while (n % p)
            {
                switch (p)
                {
                    case 4:
                        p = 2;
                        break;
                    case 2:
                        p = 3;
                        break;
                    default:
                        p += 2;
                        break;
                }
                if (p * p > n)  p = n;// no more factors
            }
            n /= p;
            _stageRadix.push_back(p);
            _stageRemainder.push_back(n);
        } while (n > 1);
    }

    /// Calculates the complex Discrete Fourier Transform.
    ///
    /// The size of the passed arrays must be passed in the constructor.
    /// The sum of the squares of the absolute values in the @c dst
    /// array will be @c N times the sum of the squares of the absolute
    /// values in the @c src array, where @c N is the size of the array.
    /// In other words, the l_2 norm of the resulting array will be
    /// @c sqrt(N) times as big as the l_2 norm of the input array.
    /// This is also the case when the inverse flag is set in the
    /// constructor. Hence when applying the same transform twice, but with
    /// the inverse flag changed the second time, then the result will
    /// be equal to the original input times @c N.
    void transform(const cpx_type * FSrc,
                   cpx_type * FDst,
                   const std::size_t stage = 0,
                   const std::size_t fstride = 1,
                   const std::size_t in_stride = 1) const
    {
        const std::size_t p = _stageRadix[stage];
        const std::size_t m = _stageRemainder[stage];
        cpx_type *const Fout_beg = FDst;
        cpx_type *const Fout_end = FDst + p * m;

        if (m == 1)
        {
            do
            {
                *FDst = *FSrc;
                FSrc += fstride * in_stride;
            } while (++FDst != 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
                transform(FSrc, FDst, stage + 1, fstride * p, in_stride);
                FSrc += fstride * in_stride;
            } while ((FDst += m) != Fout_end);
        }

        FDst = Fout_beg;

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

private:

    void kf_bfly2(cpx_type *const Fout, const size_t fstride, const std::size_t m) const
    {
        for (std::size_t k = 0; k < m; ++k)
        {
            const cpx_type t = (Fout[m + k] * _twiddles[k * fstride]) / _scale_factor;
            Fout[m + k] = Fout[k] - t;
            Fout[k] += t;
        }
    }

    void kf_bfly3(cpx_type *Fout, const std::size_t fstride, const std::size_t m) const
    {
        std::size_t k = m;
        const std::size_t m2 = 2 * m;
        const cpx_type *tw1, *tw2;
        cpx_type scratch[5];
        const cpx_type epi3 = _twiddles[fstride * m];

        tw1 = tw2 = &_twiddles[0];

        do
        {
            scratch[1] = (Fout[m] * *tw1) / _scale_factor;
            scratch[2] = (Fout[m2] * *tw2) / _scale_factor;

            scratch[3] = scratch[1] + scratch[2];
            scratch[0] = scratch[1] - scratch[2];
            tw1 += fstride;
            tw2 += fstride * 2;

            Fout[m] = Fout[0] - (scratch[3] / 2);
            scratch[0] *= epi3.imag();
            scratch[0] /= _scale_factor;

            Fout[0] += scratch[3];

            Fout[m2] = cpx_type(Fout[m].real() + scratch[0].imag(), Fout[m].imag() - scratch[0].real());

            Fout[m] += cpx_type(-scratch[0].imag(), scratch[0].real());
            ++Fout;
        } while (--k);
    }

    void kf_bfly4(cpx_type *const Fout, const std::size_t fstride, const std::size_t m) const
    {
        cpx_type scratch[7];
        const scalar_type negative_if_inverse = _inverse ? -1 : +1;

        for (std::size_t k = 0; k < m; ++k)
        {
            scratch[0] = (Fout[k + m] * _twiddles[k * fstride]) / _scale_factor;
            scratch[1] = (Fout[k + 2 * m] * _twiddles[k * fstride * 2]) / _scale_factor;
            scratch[2] = (Fout[k + 3 * m] * _twiddles[k * fstride * 3]) / _scale_factor;
            scratch[5] = Fout[k] - scratch[1];

            Fout[k] += scratch[1];
            scratch[3] = scratch[0] + scratch[2];
            scratch[4] = scratch[0] - scratch[2];
            scratch[4] = cpx_type(scratch[4].imag() * negative_if_inverse,
                                  -scratch[4].real() * negative_if_inverse);

            Fout[k + 2 * m] = Fout[k] - scratch[3];
            Fout[k] += scratch[3];
            Fout[k + m] = scratch[5] + scratch[4];
            Fout[k + 3 * m] = scratch[5] - scratch[4];
        }
    }

    void kf_bfly5(cpx_type *const Fout, const std::size_t fstride, const std::size_t m) const
    {
        cpx_type *Fout0, *Fout1, *Fout2, *Fout3, *Fout4;
        cpx_type scratch[13];
        const cpx_type ya = _twiddles[fstride * m];
        const cpx_type yb = _twiddles[fstride * 2 * m];

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

        for (std::size_t u = 0; u < m; ++u)
        {
            scratch[0] = *Fout0;

            scratch[1] = (*Fout1 * _twiddles[u * fstride]) / _scale_factor;
            scratch[2] = (*Fout2 * _twiddles[2 * u * fstride]) / _scale_factor;
            scratch[3] = (*Fout3 * _twiddles[3 * u * fstride]) / _scale_factor;
            scratch[4] = (*Fout4 * _twiddles[4 * u * fstride]) / _scale_factor;

            scratch[7] = scratch[1] + scratch[4];
            scratch[10] = scratch[1] - scratch[4];
            scratch[8] = scratch[2] + scratch[3];
            scratch[9] = scratch[2] - scratch[3];

            *Fout0 += scratch[7];
            *Fout0 += scratch[8];

            scratch[5] = scratch[0] + (cpx_type(
                    scratch[7].real() * ya.real() + scratch[8].real() * yb.real(),
                    scratch[7].imag() * ya.real() + scratch[8].imag() * yb.real() ) / _scale_factor);

            scratch[6] = cpx_type(
                    scratch[10].imag() * ya.imag() + scratch[9].imag() * yb.imag(),
                    -scratch[10].real() * ya.imag() - scratch[9].real() * yb.imag() ) / _scale_factor;

            *Fout1 = scratch[5] - scratch[6];
            *Fout4 = scratch[5] + scratch[6];

            scratch[11] = scratch[0] + (cpx_type(
                    scratch[7].real() * yb.real() + scratch[8].real() * ya.real(),
                    scratch[7].imag() * yb.real() + scratch[8].imag() * ya.real() ) / _scale_factor);

            scratch[12] = cpx_type(
                    -scratch[10].imag() * yb.imag() + scratch[9].imag() * ya.imag(),
                    scratch[10].real() * yb.imag() - scratch[9].real() * ya.imag() ) / _scale_factor;

            *Fout2 = scratch[11] + scratch[12];
            *Fout3 = scratch[11] - scratch[12];

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

    /* perform the butterfly for one stage of a mixed radix FFT */
    void kf_bfly_generic(cpx_type * const Fout, const size_t fstride, const std::size_t m, const std::size_t p) const
    {
        const cpx_type *twiddles = &_twiddles[0];
        cpx_type scratchbuf[p];

        for (std::size_t u = 0; u < m; ++u)
        {
            std::size_t k = u;
            for (std::size_t q1 = 0; q1 < p; ++q1)
            {
                scratchbuf[q1] = Fout[k];
                k += m;
            }

            k = u;
            for (std::size_t q1 = 0; q1 < p; ++q1)
            {
                std::size_t twidx = 0;
                Fout[k] = scratchbuf[0];
                for (std::size_t q = 1; q < p; ++q)
                {
                    twidx += fstride * k;
                    if (twidx >= _nfft)
                        twidx -= _nfft;
                    Fout[k] += (scratchbuf[q] * twiddles[twidx]) / _scale_factor;
                }
                k += m;
            }
        }
    }
};

#endif


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

/* kiss_fft.h
   defines kiss_fft_scalar as either short or a float type
   and defines
   typedef struct { kiss_fft_scalar r; kiss_fft_scalar i; }kiss_fft_cpx; */

#ifndef _kiss_fft_guts_h
#define _kiss_fft_guts_h

#include "kiss_fft.h"
#include "kiss_fft_log.h"
#include <limits.h>

#define MAXFACTORS 32
/* e.g. an fft of length 128 has 4 factors
 as far as kissfft is concerned
 4*4*4*2
 */

struct kiss_fft_state{
    int nfft;
    int inverse;
    int factors[2*MAXFACTORS];
    kiss_fft_cpx twiddles[1];
};

/*
  Explanation of macros dealing with complex math:

   C_MUL(m,a,b)         : m = a*b
   C_FIXDIV( c , div )  : if a fixed point impl., c /= div. noop otherwise
   C_SUB( res, a,b)     : res = a - b
   C_SUBFROM( res , a)  : res -= a
   C_ADDTO( res , a)    : res += a
 * */
#ifdef FIXED_POINT
#include <stdint.h>
#if (FIXED_POINT==32)
# define FRACBITS 31
# define SAMPPROD int64_t
#define SAMP_MAX INT32_MAX
#define SAMP_MIN INT32_MIN
#else
# define FRACBITS 15
# define SAMPPROD int32_t
#define SAMP_MAX INT16_MAX
#define SAMP_MIN INT16_MIN
#endif

#if defined(CHECK_OVERFLOW)
#  define CHECK_OVERFLOW_OP(a,op,b)  \
    if ( (SAMPPROD)(a) op (SAMPPROD)(b) > SAMP_MAX || (SAMPPROD)(a) op (SAMPPROD)(b) < SAMP_MIN ) { \
        KISS_FFT_WARNING("overflow (%d " #op" %d) = %ld", (a),(b),(SAMPPROD)(a) op (SAMPPROD)(b)); }
#endif

#   define smul(a,b) ( (SAMPPROD)(a)*(b) )
#   define sround( x )  (kiss_fft_scalar)( ( (x) + (1<<(FRACBITS-1)) ) >> FRACBITS )

#   define S_MUL(a,b) sround( smul(a,b) )

#   define C_MUL(m,a,b) \
      do{ (m).r = sround( smul((a).r,(b).r) - smul((a).i,(b).i) ); \
          (m).i = sround( smul((a).r,(b).i) + smul((a).i,(b).r) ); }while(0)

#   define DIVSCALAR(x,k) \
    (x) = sround( smul(  x, SAMP_MAX/k ) )

#   define C_FIXDIV(c,div) \
    do {    DIVSCALAR( (c).r , div);  \
        DIVSCALAR( (c).i  , div); }while (0)

#   define C_MULBYSCALAR( c, s ) \
    do{ (c).r =  sround( smul( (c).r , s ) ) ;\
        (c).i =  sround( smul( (c).i , s ) ) ; }while(0)

#else  /* not FIXED_POINT*/

#   define S_MUL(a,b) ( (a)*(b) )
#define C_MUL(m,a,b) \
    do{ (m).r = (a).r*(b).r - (a).i*(b).i;\
        (m).i = (a).r*(b).i + (a).i*(b).r; }while(0)
#   define C_FIXDIV(c,div) /* NOOP */
#   define C_MULBYSCALAR( c, s ) \
    do{ (c).r *= (s);\
        (c).i *= (s); }while(0)
#endif

#ifndef CHECK_OVERFLOW_OP
#  define CHECK_OVERFLOW_OP(a,op,b) /* noop */
#endif

#define  C_ADD( res, a,b)\
    do { \
        CHECK_OVERFLOW_OP((a).r,+,(b).r)\
        CHECK_OVERFLOW_OP((a).i,+,(b).i)\
        (res).r=(a).r+(b).r;  (res).i=(a).i+(b).i; \
    }while(0)
#define  C_SUB( res, a,b)\
    do { \
        CHECK_OVERFLOW_OP((a).r,-,(b).r)\
        CHECK_OVERFLOW_OP((a).i,-,(b).i)\
        (res).r=(a).r-(b).r;  (res).i=(a).i-(b).i; \
    }while(0)
#define C_ADDTO( res , a)\
    do { \
        CHECK_OVERFLOW_OP((res).r,+,(a).r)\
        CHECK_OVERFLOW_OP((res).i,+,(a).i)\
        (res).r += (a).r;  (res).i += (a).i;\
    }while(0)

#define C_SUBFROM( res , a)\
    do {\
        CHECK_OVERFLOW_OP((res).r,-,(a).r)\
        CHECK_OVERFLOW_OP((res).i,-,(a).i)\
        (res).r -= (a).r;  (res).i -= (a).i; \
    }while(0)

#ifdef FIXED_POINT
#  define KISS_FFT_COS(phase)  floor(.5+SAMP_MAX * cos (phase))
#  define KISS_FFT_SIN(phase)  floor(.5+SAMP_MAX * sin (phase))
#  define HALF_OF(x) ((x)>>1)
#elif defined(USE_SIMD)
#  define KISS_FFT_COS(phase) _mm_set1_ps( cos(phase) )
#  define KISS_FFT_SIN(phase) _mm_set1_ps( sin(phase) )
#  define HALF_OF(x) ((x)*_mm_set1_ps(.5))
#else
#  define KISS_FFT_COS(phase) (kiss_fft_scalar) cos(phase)
#  define KISS_FFT_SIN(phase) (kiss_fft_scalar) sin(phase)
#  define HALF_OF(x) ((x)*((kiss_fft_scalar).5))
#endif

#define  kf_cexp(x,phase) \
    do{ \
        (x)->r = KISS_FFT_COS(phase);\
        (x)->i = KISS_FFT_SIN(phase);\
    }while(0)

/* a debugging function */
#define pcpx(c)\
    KISS_FFT_DEBUG("%g + %gi\n",(double)((c)->r),(double)((c)->i))

#ifdef KISS_FFT_USE_ALLOCA
// define this to allow use of alloca instead of malloc for temporary buffers
// Temporary buffers are used in two case:
// 1. FFT sizes that have "bad" factors. i.e. not 2,3 and 5
// 2. "in-place" FFTs.  Notice the quotes, since kissfft does not really do an in-place transform.
#include <alloca.h>
#define  KISS_FFT_TMP_ALLOC(nbytes) alloca(nbytes)
#define  KISS_FFT_TMP_FREE(ptr)
#else
#define  KISS_FFT_TMP_ALLOC(nbytes) KISS_FFT_MALLOC(nbytes)
#define  KISS_FFT_TMP_FREE(ptr) KISS_FFT_FREE(ptr)
#endif

#endif /* _kiss_fft_guts_h */


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

#ifndef kiss_fft_log_h
#define kiss_fft_log_h

#define ERROR 1
#define WARNING 2
#define INFO 3
#define DEBUG 4

#define STRINGIFY(x) #x
#define TOSTRING(x) STRINGIFY(x)

#if defined(NDEBUG)
# define KISS_FFT_LOG_MSG(severity, ...) ((void)0)
#else
# define KISS_FFT_LOG_MSG(severity, ...) \
    fprintf(stderr, "[" #severity "] " __FILE__ ":" TOSTRING(__LINE__) " "); \
    fprintf(stderr, __VA_ARGS__); \
    fprintf(stderr, "\n")
#endif

#define KISS_FFT_ERROR(...) KISS_FFT_LOG_MSG(ERROR, __VA_ARGS__)
#define KISS_FFT_WARNING(...) KISS_FFT_LOG_MSG(WARNING, __VA_ARGS__)
#define KISS_FFT_INFO(...) KISS_FFT_LOG_MSG(INFO, __VA_ARGS__)
#define KISS_FFT_DEBUG(...) KISS_FFT_LOG_MSG(DEBUG, __VA_ARGS__)

#endif /* kiss_fft_log_h */


text:
cmake_minimum_required(VERSION 2.8)
project(kissfft)

add_library(kissfft
  kiss_fft.c)

target_include_directories(kissfft PUBLIC
    $<BUILD_INTERFACE:${CMAKE_CURRENT_SOURCE_DIR}>
    $<INSTALL_INTERFACE:.>)


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

#ifndef KISS_FFT_H
#define KISS_FFT_H

#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <string.h>

#ifdef __cplusplus
extern "C" {
#endif

/*
 ATTENTION!
 If you would like a :
 -- a utility that will handle the caching of fft objects
 -- real-only (no imaginary time component ) FFT
 -- a multi-dimensional FFT
 -- a command-line utility to perform ffts
 -- a command-line utility to perform fast-convolution filtering

 Then see kfc.h kiss_fftr.h kiss_fftnd.h fftutil.c kiss_fastfir.c
  in the tools/ directory.
*/

/* User may override KISS_FFT_MALLOC and/or KISS_FFT_FREE. */
#ifdef USE_SIMD
# include <xmmintrin.h>
# define kiss_fft_scalar __m128
# ifndef KISS_FFT_MALLOC
#  define KISS_FFT_MALLOC(nbytes) _mm_malloc(nbytes,16)
#  define KISS_FFT_ALIGN_CHECK(ptr)
#  define KISS_FFT_ALIGN_SIZE_UP(size) ((size + 15UL) & ~0xFUL)
# endif
# ifndef KISS_FFT_FREE
#  define KISS_FFT_FREE _mm_free
# endif
#else
# define KISS_FFT_ALIGN_CHECK(ptr)
# define KISS_FFT_ALIGN_SIZE_UP(size) (size)
# ifndef KISS_FFT_MALLOC
#  define KISS_FFT_MALLOC malloc
# endif
# ifndef KISS_FFT_FREE
#  define KISS_FFT_FREE free
# endif
#endif

#ifdef FIXED_POINT
#include <stdint.h>
# if (FIXED_POINT == 32)
#  define kiss_fft_scalar int32_t
# else
#  define kiss_fft_scalar int16_t
# endif
#else
# ifndef kiss_fft_scalar
/*  default is float */
#   define kiss_fft_scalar float
# endif
#endif

typedef struct {
    kiss_fft_scalar r;
    kiss_fft_scalar i;
}kiss_fft_cpx;

typedef struct kiss_fft_state* kiss_fft_cfg;

/*
 *  kiss_fft_alloc
 *
 *  Initialize a FFT (or IFFT) algorithm's cfg/state buffer.
 *
 *  typical usage:      kiss_fft_cfg mycfg=kiss_fft_alloc(1024,0,NULL,NULL);
 *
 *  The return value from fft_alloc is a cfg buffer used internally
 *  by the fft routine or NULL.
 *
 *  If lenmem is NULL, then kiss_fft_alloc will allocate a cfg buffer using malloc.
 *  The returned value should be free()d when done to avoid memory leaks.
 *
 *  The state can be placed in a user supplied buffer 'mem':
 *  If lenmem is not NULL and mem is not NULL and *lenmem is large enough,
 *      then the function places the cfg in mem and the size used in *lenmem
 *      and returns mem.
 *
 *  If lenmem is not NULL and ( mem is NULL or *lenmem is not large enough),
 *      then the function returns NULL and places the minimum cfg
 *      buffer size in *lenmem.
 * */

kiss_fft_cfg kiss_fft_alloc(int nfft,int inverse_fft,void * mem,size_t * lenmem);

/*
 * kiss_fft(cfg,in_out_buf)
 *
 * Perform an FFT on a complex input buffer.
 * for a forward FFT,
 * fin should be  f[0] , f[1] , ... ,f[nfft-1]
 * fout will be   F[0] , F[1] , ... ,F[nfft-1]
 * Note that each element is complex and can be accessed like
    f[k].r and f[k].i
 * */
void kiss_fft(kiss_fft_cfg cfg,const kiss_fft_cpx *fin,kiss_fft_cpx *fout);

/*
 A more generic version of the above function. It reads its input from every Nth sample.
 * */
void kiss_fft_stride(kiss_fft_cfg cfg,const kiss_fft_cpx *fin,kiss_fft_cpx *fout,int fin_stride);

/* If kiss_fft_alloc allocated a buffer, it is one contiguous
   buffer and can be simply free()d when no longer needed*/
#define kiss_fft_free KISS_FFT_FREE

/*
 Cleans up some memory that gets managed internally. Not necessary to call, but it might clean up
 your compiler output to call this before you exit.
*/
void kiss_fft_cleanup(void);

/*
 * Returns the smallest integer k, such that k>=n and k has only "fast" factors (2,3,5)
 */
int kiss_fft_next_fast_size(int n);

/* for real ffts, we need an even size */
#define kiss_fftr_next_fast_size_real(n) \
        (kiss_fft_next_fast_size( ((n)+1)>>1)<<1)

#ifdef __cplusplus
}
#endif

#endif


text:
1.3.0 2012-07-18
  removed non-standard malloc.h from kiss_fft.h

  moved -lm to end of link line

  checked various return values

  converted python Numeric code to NumPy

  fixed test of int32_t on 64 bit OS

  added padding in a couple of places to allow SIMD alignment of structs

1.2.9 2010-05-27
  threadsafe ( including OpenMP )

  first edition of kissfft.hh the C++ template fft engine

1.2.8
  Changed memory.h to string.h -- apparently more standard

  Added openmp extensions.  This can have fairly linear speedups for larger FFT sizes.

1.2.7
  Shrank the real-fft memory footprint. Thanks to Galen Seitz.

1.2.6 (Nov 14, 2006) The "thanks to GenArts" release.
  Added multi-dimensional real-optimized FFT, see tools/kiss_fftndr
  Thanks go to GenArts, Inc. for sponsoring the development.

1.2.5 (June 27, 2006) The "release for no good reason" release.
   Changed some harmless code to make some compilers' warnings go away.
   Added some more digits to pi -- why not.
   Added kiss_fft_next_fast_size() function to help people decide how much to pad.
   Changed multidimensional test from 8 dimensions to only 3 to avoid testing
   problems with fixed point (sorry Buckaroo Banzai).

1.2.4 (Oct 27, 2005)   The "oops, inverse fixed point real fft was borked" release.
   Fixed scaling bug for inverse fixed point real fft -- also fixed test code that should've been failing.
    Thanks to Jean-Marc Valin for bug report.

   Use sys/types.h for more portable types than short,int,long => int16_t,int32_t,int64_t
   If your system does not have these, you may need to define them -- but at least it breaks in a
   loud and easily fixable way -- unlike silently using the wrong size type.

   Hopefully tools/psdpng.c is fixed -- thanks to Steve Kellog for pointing out the weirdness.

1.2.3 (June 25, 2005)   The "you want to use WHAT as a sample" release.
    Added ability to use 32 bit fixed point samples -- requires a 64 bit intermediate result, a la 'long long'

    Added ability to do 4 FFTs in parallel by using SSE SIMD instructions. This is accomplished by
    using the __m128 (vector of 4 floats) as kiss_fft_scalar.  Define USE_SIMD to use this.

    I know, I know ...  this is drifting a bit from the "kiss" principle, but the speed advantages
    make it worth it for some.  Also recent gcc makes it SOO easy to use vectors of 4 floats like a POD type.

1.2.2 (May 6, 2005)   The Matthew release
    Replaced fixed point division with multiply&shift.  Thanks to Jean-Marc Valin for
    discussions regarding.  Considerable speedup for fixed-point.

    Corrected overflow protection in real fft routines  when using fixed point.
    Finder's Credit goes to Robert Oschler of robodance for pointing me at the bug.
    This also led to the CHECK_OVERFLOW_OP macro.

1.2.1 (April 4, 2004)
    compiles cleanly with just about every -W warning flag under the sun

    reorganized kiss_fft_state so it could be read-only/const. This may be useful for embedded systems
    that are willing to predeclare twiddle factors, factorization.

    Fixed C_MUL,S_MUL on 16-bit platforms.

    tmpbuf will only be allocated if input & output buffers are same
    scratchbuf will only be allocated for ffts that are not multiples of 2,3,5

    NOTE: The tmpbuf,scratchbuf changes may require synchronization code for multi-threaded apps.

1.2 (Feb 23, 2004)
    interface change -- cfg object is forward declaration of struct instead of void*
    This maintains type saftey and lets the compiler warn/error about stupid mistakes.
            (prompted by suggestion from Erik de Castro Lopo)

    small speed improvements

    added psdpng.c -- sample utility that will create png spectrum "waterfalls" from an input file
        ( not terribly useful yet)

1.1.1 (Feb 1, 2004 )
    minor bug fix -- only affects odd rank, in-place, multi-dimensional FFTs

1.1 : (Jan 30,2004)
    split sample_code/ into test/ and tools/

    Removed 2-D fft and added N-D fft (arbitrary)

    modified fftutil.c to allow multi-d FFTs

    Modified core fft routine to allow an input stride via kiss_fft_stride()
    (eased support of multi-D ffts)

    Added fast convolution filtering (FIR filtering using overlap-scrap method, with tail scrap)

    Add kfc.[ch]: the KISS FFT Cache. It takes care of allocs for you ( suggested by Oscar Lesta ).

1.0.1 (Dec 15, 2003)
    fixed bug that occurred when nfft==1. Thanks to Steven Johnson.

1.0 : (Dec 14, 2003)
    changed kiss_fft function from using a single buffer, to two buffers.
    If the same buffer pointer is supplied for both in and out, kiss will
    manage the buffer copies.

    added kiss_fft2d and kiss_fftr as separate source files (declarations in kiss_fft.h )

0.4 :(Nov 4,2003) optimized for radix 2,3,4,5

0.3 :(Oct 28, 2003) woops, version 2 didn't actually factor out any radices other than 2.
        Thanks to Steven Johnson for finding this one.

0.2 :(Oct 27, 2003) added mixed radix, only radix 2,4 optimized versions

0.1 :(May 19 2003)  initial release, radix 2 only


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

        #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:
/*
 *  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.
 */

#ifndef KISSFFT_CLASS_HH
#define KISSFFT_CLASS_HH
#include <complex>
#include <utility>
#include <vector>

template <typename scalar_t>
class kissfft
{
    public:

        typedef std::complex<scalar_t> cpx_t;

        kissfft( const std::size_t nfft,
                 const bool inverse )
            :_nfft(nfft)
            ,_inverse(inverse)
        {
            // fill twiddle factors
            _twiddles.resize(_nfft);
            const scalar_t phinc =  (_inverse?2:-2)* std::acos( (scalar_t) -1)  / _nfft;
            for (std::size_t i=0;i<_nfft;++i)
                _twiddles[i] = std::exp( cpx_t(0,i*phinc) );

            //factorize
            //start factoring out 4's, then 2's, then 3,5,7,9,...
            std::size_t n= _nfft;
            std::size_t p=4;
            do {
                while (n % p) {
                    switch (p) {
                        case 4: p = 2; break;
                        case 2: p = 3; break;
                        default: p += 2; break;
                    }
                    if (p*p>n)
                        p = n;// no more factors
                }
                n /= p;
                _stageRadix.push_back(p);
                _stageRemainder.push_back(n);
            }while(n>1);
        }

        /// Changes the FFT-length and/or the transform direction.
        ///
        /// @post The @c kissfft object will be in the same state as if it
        /// had been newly constructed with the passed arguments.
        /// However, the implementation may be faster than constructing a
        /// new fft object.
        void assign( const std::size_t nfft,
                     const bool inverse )
        {
            if ( nfft != _nfft )
            {
                kissfft tmp( nfft, inverse ); // O(n) time.
                std::swap( tmp, *this ); // this is O(1) in C++11, O(n) otherwise.
            }
            else if ( inverse != _inverse )
            {
                // conjugate the twiddle factors.
                for ( typename std::vector<cpx_t>::iterator it = _twiddles.begin();
                      it != _twiddles.end(); ++it )
                    it->imag( -it->imag() );
            }
        }

        /// Calculates the complex Discrete Fourier Transform.
        ///
        /// The size of the passed arrays must be passed in the constructor.
        /// The sum of the squares of the absolute values in the @c dst
        /// array will be @c N times the sum of the squares of the absolute
        /// values in the @c src array, where @c N is the size of the array.
        /// In other words, the l_2 norm of the resulting array will be
        /// @c sqrt(N) times as big as the l_2 norm of the input array.
        /// This is also the case when the inverse flag is set in the
        /// constructor. Hence when applying the same transform twice, but with
        /// the inverse flag changed the second time, then the result will
        /// be equal to the original input times @c N.
        void transform(const cpx_t * fft_in, cpx_t * fft_out, const std::size_t stage = 0, const std::size_t fstride = 1, const std::size_t in_stride = 1) const
        {
            const std::size_t p = _stageRadix[stage];
            const std::size_t m = _stageRemainder[stage];
            cpx_t * const Fout_beg = fft_out;
            cpx_t * const Fout_end = fft_out + p*m;

            if (m==1) {
                do{
                    *fft_out = *fft_in;
                    fft_in += fstride*in_stride;
                }while(++fft_out != 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
                    transform(fft_in, fft_out, stage+1, fstride*p,in_stride);
                    fft_in += fstride*in_stride;
                }while( (fft_out += m) != Fout_end );
            }

            fft_out=Fout_beg;

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

        /// Calculates the Discrete Fourier Transform (DFT) of a real input
        /// of size @c 2*N.
        ///
        /// The 0-th and N-th value of the DFT are real numbers. These are
        /// stored in @c dst[0].real() and @c dst[1].imag() respectively.
        /// The remaining DFT values up to the index N-1 are stored in
        /// @c dst[1] to @c dst[N-1].
        /// The other half of the DFT values can be calculated from the
        /// symmetry relation
        ///     @code
        ///         DFT(src)[2*N-k] == conj( DFT(src)[k] );
        ///     @endcode
        /// The same scaling factors as in @c transform() apply.
        ///
        /// @note For this to work, the types @c scalar_t and @c cpx_t
        /// must fulfill the following requirements:
        ///
        /// For any object @c z of type @c cpx_t,
        /// @c reinterpret_cast<scalar_t(&)[2]>(z)[0] is the real part of @c z and
        /// @c reinterpret_cast<scalar_t(&)[2]>(z)[1] is the imaginary part of @c z.
        /// For any pointer to an element of an array of @c cpx_t named @c p
        /// and any valid array index @c i, @c reinterpret_cast<T*>(p)[2*i]
        /// is the real part of the complex number @c p[i], and
        /// @c reinterpret_cast<T*>(p)[2*i+1] is the imaginary part of the
        /// complex number @c p[i].
        ///
        /// Since C++11, these requirements are guaranteed to be satisfied for
        /// @c scalar_ts being @c float, @c double or @c long @c double
        /// together with @c cpx_t being @c std::complex<scalar_t>.
        void transform_real( const scalar_t * const src,
                             cpx_t * const dst ) const
        {
            const std::size_t N = _nfft;
            if ( N == 0 )
                return;

            // perform complex FFT
            transform( reinterpret_cast<const cpx_t*>(src), dst );

            // post processing for k = 0 and k = N
            dst[0] = cpx_t( dst[0].real() + dst[0].imag(),
                               dst[0].real() - dst[0].imag() );

            // post processing for all the other k = 1, 2, ..., N-1
            const scalar_t pi = std::acos( (scalar_t) -1);
            const scalar_t half_phi_inc = ( _inverse ? pi : -pi ) / N;
            const cpx_t twiddle_mul = std::exp( cpx_t(0, half_phi_inc) );
            for ( std::size_t k = 1; 2*k < N; ++k )
            {
                const cpx_t w = (scalar_t)0.5 * cpx_t(
                     dst[k].real() + dst[N-k].real(),
                     dst[k].imag() - dst[N-k].imag() );
                const cpx_t z = (scalar_t)0.5 * cpx_t(
                     dst[k].imag() + dst[N-k].imag(),
                    -dst[k].real() + dst[N-k].real() );
                const cpx_t twiddle =
                    k % 2 == 0 ?
                    _twiddles[k/2] :
                    _twiddles[k/2] * twiddle_mul;
                dst[  k] =       w + twiddle * z;
                dst[N-k] = std::conj( w - twiddle * z );
            }
            if ( N % 2 == 0 )
                dst[N/2] = std::conj( dst[N/2] );
        }

    private:

        void kf_bfly2( cpx_t * Fout, const size_t fstride, const std::size_t m) const
        {
            for (std::size_t k=0;k<m;++k) {
                const cpx_t t = Fout[m+k] * _twiddles[k*fstride];
                Fout[m+k] = Fout[k] - t;
                Fout[k] += t;
            }
        }

        void kf_bfly3( cpx_t * Fout, const std::size_t fstride, const std::size_t m) const
        {
            std::size_t k=m;
            const std::size_t m2 = 2*m;
            const cpx_t *tw1,*tw2;
            cpx_t scratch[5];
            const cpx_t epi3 = _twiddles[fstride*m];

            tw1=tw2=&_twiddles[0];

            do{
                scratch[1] = Fout[m]  * *tw1;
                scratch[2] = Fout[m2] * *tw2;

                scratch[3] = scratch[1] + scratch[2];
                scratch[0] = scratch[1] - scratch[2];
                tw1 += fstride;
                tw2 += fstride*2;

                Fout[m] = Fout[0] - scratch[3]*scalar_t(0.5);
                scratch[0] *= epi3.imag();

                Fout[0] += scratch[3];

                Fout[m2] = cpx_t(  Fout[m].real() + scratch[0].imag() , Fout[m].imag() - scratch[0].real() );

                Fout[m] += cpx_t( -scratch[0].imag(),scratch[0].real() );
                ++Fout;
            }while(--k);
        }

        void kf_bfly4( cpx_t * const Fout, const std::size_t fstride, const std::size_t m) const
        {
            cpx_t scratch[7];
            const scalar_t negative_if_inverse = _inverse ? -1 : +1;
            for (std::size_t k=0;k<m;++k) {
                scratch[0] = Fout[k+  m] * _twiddles[k*fstride  ];
                scratch[1] = Fout[k+2*m] * _twiddles[k*fstride*2];
                scratch[2] = Fout[k+3*m] * _twiddles[k*fstride*3];
                scratch[5] = Fout[k] - scratch[1];

                Fout[k] += scratch[1];
                scratch[3] = scratch[0] + scratch[2];
                scratch[4] = scratch[0] - scratch[2];
                scratch[4] = cpx_t( scratch[4].imag()*negative_if_inverse ,
                                      -scratch[4].real()*negative_if_inverse );

                Fout[k+2*m]  = Fout[k] - scratch[3];
                Fout[k    ]+= scratch[3];
                Fout[k+  m] = scratch[5] + scratch[4];
                Fout[k+3*m] = scratch[5] - scratch[4];
            }
        }

        void kf_bfly5( cpx_t * const Fout, const std::size_t fstride, const std::size_t m) const
        {
            cpx_t *Fout0,*Fout1,*Fout2,*Fout3,*Fout4;
            cpx_t scratch[13];
            const cpx_t ya = _twiddles[fstride*m];
            const cpx_t yb = _twiddles[fstride*2*m];

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

            for ( std::size_t u=0; u<m; ++u ) {
                scratch[0] = *Fout0;

                scratch[1] = *Fout1 * _twiddles[  u*fstride];
                scratch[2] = *Fout2 * _twiddles[2*u*fstride];
                scratch[3] = *Fout3 * _twiddles[3*u*fstride];
                scratch[4] = *Fout4 * _twiddles[4*u*fstride];

                scratch[7] = scratch[1] + scratch[4];
                scratch[10]= scratch[1] - scratch[4];
                scratch[8] = scratch[2] + scratch[3];
                scratch[9] = scratch[2] - scratch[3];

                *Fout0 += scratch[7];
                *Fout0 += scratch[8];

                scratch[5] = scratch[0] + cpx_t(
                        scratch[7].real()*ya.real() + scratch[8].real()*yb.real(),
                        scratch[7].imag()*ya.real() + scratch[8].imag()*yb.real()
                        );

                scratch[6] =  cpx_t(
                         scratch[10].imag()*ya.imag() + scratch[9].imag()*yb.imag(),
                        -scratch[10].real()*ya.imag() - scratch[9].real()*yb.imag()
                        );

                *Fout1 = scratch[5] - scratch[6];
                *Fout4 = scratch[5] + scratch[6];

                scratch[11] = scratch[0] +
                    cpx_t(
                            scratch[7].real()*yb.real() + scratch[8].real()*ya.real(),
                            scratch[7].imag()*yb.real() + scratch[8].imag()*ya.real()
                            );

                scratch[12] = cpx_t(
                        -scratch[10].imag()*yb.imag() + scratch[9].imag()*ya.imag(),
                         scratch[10].real()*yb.imag() - scratch[9].real()*ya.imag()
                        );

                *Fout2 = scratch[11] + scratch[12];
                *Fout3 = scratch[11] - scratch[12];

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

        /* perform the butterfly for one stage of a mixed radix FFT */
        void kf_bfly_generic(
                cpx_t * const Fout,
                const size_t fstride,
                const std::size_t m,
                const std::size_t p
                ) const
        {
            const cpx_t * twiddles = &_twiddles[0];

            if(p > _scratchbuf.size()) _scratchbuf.resize(p);

            for ( std::size_t u=0; u<m; ++u ) {
                std::size_t k = u;
                for ( std::size_t q1=0 ; q1<p ; ++q1 ) {
                    _scratchbuf[q1] = Fout[ k  ];
                    k += m;
                }

                k=u;
                for ( std::size_t q1=0 ; q1<p ; ++q1 ) {
                    std::size_t twidx=0;
                    Fout[ k ] = _scratchbuf[0];
                    for ( std::size_t q=1;q<p;++q ) {
                        twidx += fstride * k;
                        if (twidx>=_nfft)
                          twidx-=_nfft;
                        Fout[ k ] += _scratchbuf[q] * twiddles[twidx];
                    }
                    k += m;
                }
            }
        }

        std::size_t _nfft;
        bool _inverse;
        std::vector<cpx_t> _twiddles;
        std::vector<std::size_t> _stageRadix;
        std::vector<std::size_t> _stageRemainder;
        mutable std::vector<cpx_t> _scratchbuf;
};
#endif


text:
If you are reading this, it means you think you may be interested in using the SIMD extensions in kissfft
to do 4 *separate* FFTs at once.

Beware! Beyond here there be dragons!

This API is not easy to use, is not well documented, and breaks the KISS principle.

Still reading? Okay, you may get rewarded for your patience with a considerable speedup
(2-3x) on intel x86 machines with SSE if you are willing to jump through some hoops.

The basic idea is to use the packed 4 float __m128 data type as a scalar element.
This means that the format is pretty convoluted. It performs 4 FFTs per fft call on signals A,B,C,D.

For complex data, the data is interlaced as follows:
rA0,rB0,rC0,rD0,      iA0,iB0,iC0,iD0,   rA1,rB1,rC1,rD1, iA1,iB1,iC1,iD1 ...
where "rA0" is the real part of the zeroth sample for signal A

Real-only data is laid out:
rA0,rB0,rC0,rD0,     rA1,rB1,rC1,rD1,      ...

Compile with gcc flags something like
-O3 -mpreferred-stack-boundary=4  -DUSE_SIMD=1 -msse

Be aware of SIMD alignment.  This is the most likely cause of segfaults.
The code within kissfft uses scratch variables on the stack.
With SIMD, these must have addresses on 16 byte boundaries.
Search on "SIMD alignment" for more info.

Robin at Divide Concept was kind enough to share his code for formatting to/from the SIMD kissfft.
I have not run it -- use it at your own risk.  It appears to do 4xN and Nx4 transpositions
(out of place).

void SSETools::pack128(float* target, float* source, unsigned long size128)
{
   __m128* pDest = (__m128*)target;
   __m128* pDestEnd = pDest+size128;
   float* source0=source;
   float* source1=source0+size128;
   float* source2=source1+size128;
   float* source3=source2+size128;

   while(pDest<pDestEnd)
   {
       *pDest=_mm_set_ps(*source3,*source2,*source1,*source0);
       source0++;
       source1++;
       source2++;
       source3++;
       pDest++;
   }
}

void SSETools::unpack128(float* target, float* source, unsigned long size128)
{

   float* pSrc = source;
   float* pSrcEnd = pSrc+size128*4;
   float* target0=target;
   float* target1=target0+size128;
   float* target2=target1+size128;
   float* target3=target2+size128;

   while(pSrc<pSrcEnd)
   {
       *target0=pSrc[0];
       *target1=pSrc[1];
       *target2=pSrc[2];
       *target3=pSrc[3];
       target0++;
       target1++;
       target2++;
       target3++;
       pSrc+=4;
   }
}


text:
KFVER=131

DATATYPE ?= float

PREFIX ?= /usr/local
LIBDIR ?= $(PREFIX)/lib

INSTALL ?= install

ifeq ($(shell uname -s),Darwin)
	SHARED_NAME := libkissfft.dylib
	SHARED_FLAGS := -Wl,-install_name,$(SHARED_NAME)
else
	SHARED_NAME := libkissfft.so
	SHARED_FLAGS := -Wl,-soname,$(SHARED_NAME)
endif

all:
	gcc -Wall -fPIC -c *.c -Dkiss_fft_scalar=$(DATATYPE) -o kiss_fft.o
	ar crus libkissfft.a kiss_fft.o
	gcc -shared $(SHARED_FLAGS) -o $(SHARED_NAME) kiss_fft.o

install: all
	$(INSTALL) -Dt $(LIBDIR) $(SHARED_NAME)

doc:
	@echo "Start by reading the README file.  If you want to build and test lots of stuff, do a 'make testall'"
	@echo "but be aware that 'make testall' has dependencies that the basic kissfft software does not."
	@echo "It is generally unneeded to run these tests yourself, unless you plan on changing the inner workings"
	@echo "of kissfft and would like to make use of its regression tests."

testall:
	# The simd and int32_t types may or may not work on your machine
	make -C test testcpp && test/testcpp
	make -C test DATATYPE=simd CFLAGADD="$(CFLAGADD)" test
	make -C test DATATYPE=int32_t CFLAGADD="$(CFLAGADD)" test
	make -C test DATATYPE=int16_t CFLAGADD="$(CFLAGADD)" test
	make -C test DATATYPE=float CFLAGADD="$(CFLAGADD)" test
	make -C test DATATYPE=double CFLAGADD="$(CFLAGADD)" test
	make -C test testsse
	echo "all tests passed"

tarball: clean
	git archive --prefix=kissfft/ -o kissfft$(KFVER).tar.gz v$(KFVER)
	git archive --prefix=kissfft/ -o kissfft$(KFVER).zip v$(KFVER)

clean:
	cd test && make clean
	cd tools && make clean
	rm -f kiss_fft*.tar.gz *~ *.pyc kiss_fft*.zip

asm: kiss_fft.s

kiss_fft.s: kiss_fft.c kiss_fft.h _kiss_fft_guts.h
	[ -e kiss_fft.s ] && mv kiss_fft.s kiss_fft.s~ || true
	gcc -S kiss_fft.c -O3 -mtune=native -ffast-math -fomit-frame-pointer -unroll-loops -dA -fverbose-asm
	gcc -o kiss_fft_short.s -S kiss_fft.c -O3 -mtune=native -ffast-math -fomit-frame-pointer -dA -fverbose-asm -DFIXED_POINT
	[ -e kiss_fft.s~ ] && diff kiss_fft.s~ kiss_fft.s || true


text:
Copyright (c) 2003-2010 Mark Borgerding . All rights reserved.

KISS FFT is provided under:

  SPDX-License-Identifier: BSD-3-Clause

Being under the terms of the BSD 3-clause "New" or "Revised" License,
according with:

  LICENSES/BSD-3-Clause