CDSPRealFFT.h

Go to the documentation of this file.

//$ nobt
//$ nocpp

 
#ifndef R8B_CDSPREALFFT_INCLUDED
#define R8B_CDSPREALFFT_INCLUDED
 
#include "r8bbase.h"
 
#if R8B_PFFFT_DOUBLE
    #include "fft/pffft_double.h"
#elif R8B_PFFFT
    #include "fft/pffft.h"
#elif !R8B_IPP
    #include "fft/fft4g.h"
#endif // !R8B_IPP
 
namespace r8b {
 
#if R8B_PFFFT
    typedef float realfft_t; 
#else // R8B_PFFFT
    typedef double realfft_t; 
#endif // R8B_PFFFT

 
class CDSPRealFFTSetup : public R8B_BASECLASS
{
    R8BNOCTOR( CDSPRealFFTSetup )
 
    friend class CPtrKeeper< CDSPRealFFTSetup >;
    friend class CDSPRealFFT;
    friend class CDSPRealFFTKeeper;
 
public:
 
    double getInvMulConst() const
    {
        return( InvMulConst );
    }

 
    int getLenBits() const
    {
        return( LenBits );
    }

 
    int getLen() const
    {
        return( Len );
    }
 
private:
    int LenBits; 
    int Len; 
    double InvMulConst; 
 
    #if R8B_IPP
        IppsFFTSpec_R_64f* SpecPtr; 
        CFixedBuffer< unsigned char > SpecBuffer; 
        int WorkBufferSize; 
    #elif R8B_PFFFT
        PFFFT_Setup* setup; 
    #elif R8B_PFFFT_DOUBLE
        PFFFTD_Setup* setup; 
    #else // R8B_PFFFT_DOUBLE
        CFixedBuffer< int > wi; 
        CFixedBuffer< double > wd; 
    #endif // R8B_PFFFT_DOUBLE

 
    CDSPRealFFTSetup( const int aLenBits )
        : LenBits( aLenBits )
        , Len( 1 << aLenBits )
    #if R8B_OOURA
        , InvMulConst( 2.0 / Len )
    #else // R8B_OOURA
        , InvMulConst( 1.0 / Len )
    #endif // R8B_OOURA
    {
    #if R8B_IPP
 
        int SpecBufferSize;
        int InitBufferSize;
 
        ippsFFTGetSize_R_64f( LenBits, IPP_FFT_NODIV_BY_ANY,
            ippAlgHintFast, &SpecBufferSize, &InitBufferSize,
            &WorkBufferSize );
 
        CFixedBuffer< unsigned char > InitBuffer( InitBufferSize );
        SpecBuffer.alloc( SpecBufferSize );
 
        ippsFFTInit_R_64f( &SpecPtr, LenBits, IPP_FFT_NODIV_BY_ANY,
            ippAlgHintFast, SpecBuffer, InitBuffer );
 
    #elif R8B_PFFFT
 
        setup = pffft_new_setup( Len, PFFFT_REAL );
 
    #elif R8B_PFFFT_DOUBLE
 
        setup = pffftd_new_setup( Len, PFFFT_REAL );
 
    #else // R8B_PFFFT_DOUBLE
 
        wi.alloc( (int) ceil( 2.0 + sqrt( (double) ( Len >> 1 ))));
        wd.alloc( Len >> 1 );
 
        ooura_fft :: rdftinit( Len, wi, wd );
 
    #endif // R8B_PFFFT_DOUBLE
    }
 
    ~CDSPRealFFTSetup()
    {
        #if R8B_PFFFT
            pffft_destroy_setup( setup );
        #elif R8B_PFFFT_DOUBLE
            pffftd_destroy_setup( setup );
        #endif // R8B_PFFFT_DOUBLE
    }
};

 
class CDSPRealFFT : public R8B_BASECLASS,
    private CSinglyLinkedListItem< CDSPRealFFT >
{
    R8BNOCTOR( CDSPRealFFT )
 
    friend class CSinglyLinkedListItem< CDSPRealFFT >;
    friend class CPtrKeeper< CDSPRealFFT >;
    friend class CDSPRealFFTKeeper;
 
public:
 
    double getInvMulConst() const
    {
        return( s -> InvMulConst );
    }

 
    int getLenBits() const
    {
        return( s -> LenBits );
    }

 
    int getLen() const
    {
        return( Len );
    }

 
    realfft_t* forward( double* const p ) const
    {
    #if R8B_IPP
 
        ippsFFTFwd_RToPerm_64f( p, p, s -> SpecPtr, WorkBuffer );
 
        return( p );
 
    #elif R8B_PFFFT
 
        float* const tp = &work[ Len ];
        int i;
 
        for( i = 0; i < Len; i++ )
        {
            tp[ i ] = (float) p[ i ];
        }
 
        float* const op = construct_ptr< float >( (void*) p, (size_t) Len );
 
        pffft_transform( s -> setup, tp, op, work, PFFFT_FORWARD );
 
        return( op );
 
    #elif R8B_PFFFT_DOUBLE
 
        pffftd_transform( s -> setup, p, p, work, PFFFT_FORWARD );
 
        return( p );
 
    #else // R8B_PFFFT_DOUBLE
 
        ooura_fft :: rdft( Len, 1, p, s -> wi, s -> wd );
 
        return( p );
 
    #endif // R8B_PFFFT_DOUBLE
    }

 
    void inverse( realfft_t* const p ) const
    {
    #if R8B_IPP
 
        ippsFFTInv_PermToR_64f( p, p, s -> SpecPtr, WorkBuffer );
 
    #elif R8B_PFFFT
 
        float* const tp = &work[ Len ];
 
        pffft_transform( s -> setup, p, tp, work, PFFFT_BACKWARD );
 
        double* const op = construct_ptr< double >( (void*) p, (size_t) Len );
        int i;
 
        for( i = 0; i < Len; i++ )
        {
            op[ i ] = tp[ i ];
        }
 
    #elif R8B_PFFFT_DOUBLE
 
        pffftd_transform( s -> setup, p, p, work, PFFFT_BACKWARD );
 
    #else // R8B_PFFFT_DOUBLE
 
        ooura_fft :: rdft( Len, -1, p, s -> wi, s -> wd );
 
    #endif // R8B_PFFFT_DOUBLE
    }

 
    realfft_t* getWorkBuf() const
    {
    #if R8B_PFFFT || R8B_PFFFT_DOUBLE
 
        return( work );
 
    #else // R8B_PFFFT || R8B_PFFFT_DOUBLE
 
        return( R8B_NULL );
 
    #endif // R8B_PFFFT || R8B_PFFFT_DOUBLE
    }

 
    realfft_t* reorderForward( realfft_t* const p, realfft_t* const op ) const
    {
    #if R8B_PFFFT
 
        pffft_zreorder( s -> setup, p, op, PFFFT_FORWARD );
 
        return( op );
 
    #elif R8B_PFFFT_DOUBLE
 
        pffftd_zreorder( s -> setup, p, op, PFFFT_FORWARD );
 
        return( op );
 
    #else // R8B_PFFFT_DOUBLE
 
        (void) op;
 
        return( p );
 
    #endif // R8B_PFFFT_DOUBLE
    }

 
    realfft_t* reorderInverse( const realfft_t* const ip,
        double* const op ) const
    {
    #if R8B_PFFFT
 
        float* const op2 = construct_ptr< float >( (void*) op, Len );
 
        pffft_zreorder( s -> setup, ip, op2, PFFFT_BACKWARD );
 
        return( op2 );
 
    #elif R8B_PFFFT_DOUBLE
 
        pffftd_zreorder( s -> setup, ip, op, PFFFT_BACKWARD );
 
        return( op );
 
    #else // R8B_PFFFT_DOUBLE
 
        if( ip != op )
        {
            memcpy( op, ip, (size_t) Len * sizeof( double ));
        }
 
        return( op );
 
    #endif // R8B_PFFFT_DOUBLE
    }

 
    static void setBinNyquist( realfft_t* const p, const realfft_t v )
    {
    #if R8B_PFFFT || R8B_PFFFT_DOUBLE
 
        p[ 4 ] = v;
 
    #else // R8B_PFFFT || R8B_PFFFT_DOUBLE
 
        p[ 1 ] = v;
 
    #endif // R8B_PFFFT || R8B_PFFFT_DOUBLE
    }

 
    void multiplyBlocks( const realfft_t* const ip1,
        const realfft_t* const ip2, realfft_t* const op ) const
    {
    #if R8B_IPP
 
        ippsMulPerm_64f( (Ipp64f*) ip1, (Ipp64f*) ip2, (Ipp64f*) op, Len );
 
    #elif R8B_PFFFT
 
        pffft_zconvolve( s -> setup, ip1, ip2, op );
 
    #elif R8B_PFFFT_DOUBLE
 
        pffftd_zconvolve( s -> setup, ip1, ip2, op );
 
    #else // R8B_PFFFT_DOUBLE
 
        op[ 0 ] = ip1[ 0 ] * ip2[ 0 ];
        op[ 1 ] = ip1[ 1 ] * ip2[ 1 ];
 
        int i = 2;
 
        while( i < Len )
        {
            op[ i ] = ip1[ i ] * ip2[ i ] - ip1[ i + 1 ] * ip2[ i + 1 ];
            op[ i + 1 ] = ip1[ i ] * ip2[ i + 1 ] + ip1[ i + 1 ] * ip2[ i ];
            i += 2;
        }
 
    #endif // R8B_PFFFT_DOUBLE
    }

 
    void multiplyBlocks( const realfft_t* const ip, realfft_t* const op ) const
    {
    #if R8B_IPP
 
        ippsMulPerm_64f( (Ipp64f*) op, (Ipp64f*) ip, (Ipp64f*) op, Len );
 
    #elif R8B_PFFFT
 
        pffft_zconvolve( s -> setup, ip, op, op );
 
    #elif R8B_PFFFT_DOUBLE
 
        pffftd_zconvolve( s -> setup, ip, op, op );
 
    #else // R8B_PFFFT_DOUBLE
 
        op[ 0 ] *= ip[ 0 ];
        op[ 1 ] *= ip[ 1 ];
 
        realfft_t t;
        int i = 2;
 
        while( i < Len )
        {
            t = op[ i ] * ip[ i ] - op[ i + 1 ] * ip[ i + 1 ];
            op[ i + 1 ] = op[ i ] * ip[ i + 1 ] + op[ i + 1 ] * ip[ i ];
            op[ i ] = t;
            i += 2;
        }
 
    #endif // R8B_PFFFT_DOUBLE
    }

 
    void multiplyBlocksZP( const realfft_t* ip, realfft_t* op ) const
    {
    // SIMD implementations assume that pointers are address-aligned.
 
    #if R8B_PFFFT_DOUBLE
 
        pffftd_zconvolve_zp( s -> setup, op, ip, op );
 
    #elif R8B_PFFFT && defined( R8B_SSE2 )
 
        int c16 = Len >> 4;
 
        while( c16 != 0 )
        {
            const __m128 iv1 = _mm_load_ps( ip );
            const __m128 iv2 = _mm_load_ps( ip + 4 );
            const __m128 ov1 = _mm_load_ps( op );
            const __m128 ov2 = _mm_load_ps( op + 4 );
            _mm_store_ps( op, _mm_mul_ps( iv1, ov1 ));
            _mm_store_ps( op + 4, _mm_mul_ps( iv2, ov2 ));
 
            const __m128 iv3 = _mm_load_ps( ip + 8 );
            const __m128 iv4 = _mm_load_ps( ip + 12 );
            const __m128 ov3 = _mm_load_ps( op + 8 );
            const __m128 ov4 = _mm_load_ps( op + 12 );
            _mm_store_ps( op + 8, _mm_mul_ps( iv3, ov3 ));
            _mm_store_ps( op + 12, _mm_mul_ps( iv4, ov4 ));
 
            ip += 16;
            op += 16;
            c16--;
        }
 
        int c = Len & 15;
 
        while( c != 0 )
        {
            *op *= *ip;
            ip++;
            op++;
            c--;
        }
 
    #elif R8B_PFFFT && defined( R8B_NEON )
 
        int c16 = Len >> 4;
 
        while( c16 != 0 )
        {
            const float32x4_t iv1 = vld1q_f32( ip );
            const float32x4_t iv2 = vld1q_f32( ip + 4 );
            const float32x4_t ov1 = vld1q_f32( op );
            const float32x4_t ov2 = vld1q_f32( op + 4 );
            vst1q_f32( op, vmulq_f32( iv1, ov1 ));
            vst1q_f32( op + 4, vmulq_f32( iv2, ov2 ));
 
            const float32x4_t iv3 = vld1q_f32( ip + 8 );
            const float32x4_t iv4 = vld1q_f32( ip + 12 );
            const float32x4_t ov3 = vld1q_f32( op + 8 );
            const float32x4_t ov4 = vld1q_f32( op + 12 );
            vst1q_f32( op + 8, vmulq_f32( iv3, ov3 ));
            vst1q_f32( op + 12, vmulq_f32( iv4, ov4 ));
 
            ip += 16;
            op += 16;
            c16--;
        }
 
        int c = Len & 15;
 
        while( c != 0 )
        {
            *op *= *ip;
            ip++;
            op++;
            c--;
        }
 
    #elif !R8B_PFFFT && defined( R8B_SSE2 )
 
        int c8 = Len >> 3;
 
        while( c8 != 0 )
        {
            const __m128d iv1 = _mm_load_pd( ip );
            const __m128d iv2 = _mm_load_pd( ip + 2 );
            const __m128d ov1 = _mm_load_pd( op );
            const __m128d ov2 = _mm_load_pd( op + 2 );
            _mm_store_pd( op, _mm_mul_pd( iv1, ov1 ));
            _mm_store_pd( op + 2, _mm_mul_pd( iv2, ov2 ));
 
            const __m128d iv3 = _mm_load_pd( ip + 4 );
            const __m128d iv4 = _mm_load_pd( ip + 6 );
            const __m128d ov3 = _mm_load_pd( op + 4 );
            const __m128d ov4 = _mm_load_pd( op + 6 );
            _mm_store_pd( op + 4, _mm_mul_pd( iv3, ov3 ));
            _mm_store_pd( op + 6, _mm_mul_pd( iv4, ov4 ));
 
            ip += 8;
            op += 8;
            c8--;
        }
 
        int c = Len & 7;
 
        while( c != 0 )
        {
            *op *= *ip;
            ip++;
            op++;
            c--;
        }
 
    #elif !R8B_PFFFT && defined( R8B_NEON )
 
        int c8 = Len >> 3;
 
        while( c8 != 0 )
        {
            const float64x2_t iv1 = vld1q_f64( ip );
            const float64x2_t iv2 = vld1q_f64( ip + 2 );
            const float64x2_t ov1 = vld1q_f64( op );
            const float64x2_t ov2 = vld1q_f64( op + 2 );
            vst1q_f64( op, vmulq_f64( iv1, ov1 ));
            vst1q_f64( op + 2, vmulq_f64( iv2, ov2 ));
 
            const float64x2_t iv3 = vld1q_f64( ip + 4 );
            const float64x2_t iv4 = vld1q_f64( ip + 6 );
            const float64x2_t ov3 = vld1q_f64( op + 4 );
            const float64x2_t ov4 = vld1q_f64( op + 6 );
            vst1q_f64( op + 4, vmulq_f64( iv3, ov3 ));
            vst1q_f64( op + 6, vmulq_f64( iv4, ov4 ));
 
            ip += 8;
            op += 8;
            c8--;
        }
 
        int c = Len & 7;
 
        while( c != 0 )
        {
            *op *= *ip;
            ip++;
            op++;
            c--;
        }
 
    #else // SIMD
 
        int i;
 
        for( i = 0; i < Len; i++ )
        {
            op[ i ] *= ip[ i ];
        }
 
    #endif // SIMD
    }

 
    void convertToZP( realfft_t* const p ) const
    {
    #if R8B_PFFFT
 
        pffft_zreorder( s -> setup, p, work, PFFFT_FORWARD );
 
        int i = 2;
 
        while( i < Len )
        {
            work[ i + 1 ] = work[ i ];
            i += 2;
        }
 
        pffft_zreorder( s -> setup, work, p, PFFFT_BACKWARD );
 
    #elif R8B_PFFFT_DOUBLE
 
        pffftd_zconvert_zp( s -> setup, p, p, (double) Len );
 
    #else // R8B_PFFFT_DOUBLE
 
        int i = 2;
 
        while( i < Len )
        {
            p[ i + 1 ] = p[ i ];
            i += 2;
        }
 
    #endif // R8B_PFFFT_DOUBLE
    }
 
private:
    const CDSPRealFFTSetup* s; 
    int Len; 
 
    #if R8B_IPP
        CFixedBuffer< unsigned char > WorkBuffer; 
    #elif R8B_PFFFT
        CFixedBuffer< float > work; 
    #elif R8B_PFFFT_DOUBLE
        CFixedBuffer< double > work; 
    #endif // R8B_PFFFT_DOUBLE

 
    CDSPRealFFT( const CDSPRealFFTSetup* const s0 )
        : s( s0 )
        , Len( s -> Len )
    {
    #if R8B_IPP
 
        WorkBuffer.alloc( s -> WorkBufferSize );
 
    #elif R8B_PFFFT
 
        work.alloc( Len * 2 );
 
    #elif R8B_PFFFT_DOUBLE
 
        work.alloc( Len );
 
    #endif // R8B_PFFFT_DOUBLE
    }
};

 
class CDSPRealFFTKeeper : public R8B_BASECLASS
{
    R8BNOCTOR( CDSPRealFFTKeeper )
 
public:
    CDSPRealFFTKeeper()
        : Object( R8B_NULL )
    {
    }

 
    CDSPRealFFTKeeper( const int LenBits )
    {
        Object = acquire( LenBits );
    }
 
    ~CDSPRealFFTKeeper()
    {
        if( Object != R8B_NULL )
        {
            release( Object );
        }
    }

 
    const CDSPRealFFT* operator -> () const
    {
        R8BASSERT( Object != R8B_NULL );
 
        return( Object );
    }

 
    void init( const int LenBits )
    {
        if( Object != R8B_NULL )
        {
            if( Object -> getLenBits() == LenBits )
            {
                return;
            }
 
            release( Object );
        }
 
        Object = acquire( LenBits );
    }

 
    void reset()
    {
        if( Object != R8B_NULL )
        {
            release( Object );
            Object = R8B_NULL;
        }
    }
 
private:
    CDSPRealFFT* Object; 

 
    CDSPRealFFT* acquire( const int LenBits )
    {
        R8BASSERT( LenBits > 0 && LenBits <= 30 );
 
        R8BSYNC( getStateSync() );
 
        CDSPRealFFTSetup* s = getFFTSetupObjects()[ LenBits ];
 
        if( s == R8B_NULL )
        {
            s = new CDSPRealFFTSetup( LenBits );
            getFFTSetupObjects()[ LenBits ] = s;
        }
 
        #if R8B_OOURA
 
            if( getFFTObjects()[ LenBits ] == R8B_NULL )
            {
                getFFTObjects()[ LenBits ] = new CDSPRealFFT( s );
            }
 
            return( getFFTObjects()[ LenBits ]);
 
        #else // R8B_OOURA
 
            if( getFFTObjects()[ LenBits ] == R8B_NULL )
            {
                return( new CDSPRealFFT( s ));
            }
 
            CDSPRealFFT* const ffto = getFFTObjects()[ LenBits ].unkeep();
            getFFTObjects()[ LenBits ] = ffto -> Next;
            ffto -> Next = R8B_NULL;
 
            return( ffto );
 
        #endif // R8B_OOURA
    }

 
    void release( CDSPRealFFT* const ffto )
    {
    #if R8B_OOURA
 
        (void) ffto;
 
    #else // R8B_OOURA
 
        R8BSYNC( getStateSync() );
 
        ffto -> Next = getFFTObjects()[ ffto -> getLenBits() ].unkeep();
        getFFTObjects()[ ffto -> getLenBits() ] = ffto;
 
    #endif // R8B_OOURA
    }

 
    static CPtrKeeper< CDSPRealFFTSetup >* getFFTSetupObjects()
    {
        R8B_EXITDTOR static CPtrKeeper< CDSPRealFFTSetup > FFTSetups[ 31 ];
 
        return( FFTSetups );
    }

 
    static CPtrKeeper< CDSPRealFFT >* getFFTObjects()
    {
        R8B_EXITDTOR static CPtrKeeper< CDSPRealFFT > FFTObjects[ 31 ];
 
        return( FFTObjects );
    }

 
    static CSyncObject& getStateSync()
    {
        R8B_EXITDTOR static CSyncObject StateSync;
 
        return( StateSync );
    }
};

 
inline void calcMinPhaseTransform( double* const Kernel, const int KernelLen,
    const int LenMult = 2, const bool DoFinalMul = true,
    double* const DCGroupDelay = R8B_NULL )
{
    R8BASSERT( KernelLen > 0 );
    R8BASSERT( LenMult >= 2 );
 
    const int LenBits = getBitOccupancy(( KernelLen * LenMult ) - 1 );
    const int Len = 1 << LenBits;
    const int Len2 = Len >> 1;
    int i;
 
    CFixedBuffer< double > ip( Len );
    CFixedBuffer< realfft_t > ip2( Len2 + 1 );
 
    memcpy( &ip[ 0 ], Kernel, (size_t) KernelLen * sizeof( ip[ 0 ]));
    memset( &ip[ KernelLen ], 0,
        (size_t) ( Len - KernelLen ) * sizeof( ip[ 0 ]));
 
    CDSPRealFFTKeeper ffto( LenBits );
 
    realfft_t* aip = ffto -> reorderForward( ffto -> forward( ip ),
        ffto -> getWorkBuf() );
 
    realfft_t* const aip2 = &ip2[ 0 ];
 
    const realfft_t nzbias = ( sizeof( realfft_t ) < sizeof( double ) ?
        1e-35 : 1e-300 );
 
    // Create the "log |c|" spectrum while saving the original power spectrum
    // in the "ip2" buffer.
 
    aip2[ 0 ] = aip[ 0 ];
    aip[ 0 ] = log( fabs( aip[ 0 ]) + nzbias );
    aip2[ Len2 ] = aip[ 1 ];
    aip[ 1 ] = log( fabs( aip[ 1 ]) + nzbias );
 
    for( i = 1; i < Len2; i++ )
    {
        aip2[ i ] = sqrt( aip[ i * 2 ] * aip[ i * 2 ] +
            aip[ i * 2 + 1 ] * aip[ i * 2 + 1 ]);
 
        aip[ i * 2 ] = log( aip2[ i ] + nzbias );
        aip[ i * 2 + 1 ] = 0.0;
    }
 
    // Convert to cepstrum and apply discrete Hilbert transform.
 
    ffto -> inverse( ffto -> reorderInverse( aip, ip ));
 
    const double m1 = ffto -> getInvMulConst();
    const double m2 = -m1;
 
    ip[ 0 ] = 0.0;
 
    for( i = 1; i < Len2; i++ )
    {
        ip[ i ] *= m1;
    }
 
    ip[ Len2 ] = 0.0;
 
    for( i = Len2 + 1; i < Len; i++ )
    {
        ip[ i ] *= m2;
    }
 
    // Convert Hilbert-transformed cepstrum back to the "log |c|" spectrum and
    // perform its exponentiation, multiplied by the power spectrum previously
    // saved in the "ip2" buffer.
 
    aip = ffto -> reorderForward( ffto -> forward( ip ),
        ffto -> getWorkBuf() );
 
    aip[ 0 ] = aip2[ 0 ];
    aip[ 1 ] = aip2[ Len2 ];
 
    for( i = 1; i < Len2; i++ )
    {
        aip[ i * 2 + 0 ] = cos( aip[ i * 2 + 1 ]) * aip2[ i ];
        aip[ i * 2 + 1 ] = sin( aip[ i * 2 + 1 ]) * aip2[ i ];
    }
 
    ffto -> inverse( ffto -> reorderInverse( aip, ip ));
 
    if( DoFinalMul )
    {
        for( i = 0; i < KernelLen; i++ )
        {
            Kernel[ i ] = ip[ i ] * m1;
        }
    }
    else
    {
        memcpy( &Kernel[ 0 ], &ip[ 0 ],
            (size_t) KernelLen * sizeof( Kernel[ 0 ]));
    }
 
    if( DCGroupDelay != R8B_NULL )
    {
        *DCGroupDelay = calcFIRFilterGroupDelay( Kernel, KernelLen, 0.0 );
    }
}
 
} // namespace r8b
 
#endif // R8B_CDSPREALFFT_INCLUDED