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 "pffft_double/pffft_double.h"
#elif R8B_PFFFT
    #include "pffft.h"
#elif !R8B_IPP
    #include "fft4g.h"
#endif // !R8B_IPP
 
namespace r8b {
 
class CDSPRealFFT : public R8B_BASECLASS
{
    R8BNOCTOR( CDSPRealFFT );
 
    friend class CDSPRealFFTKeeper;
 
public:
    double getInvMulConst() const
    {
        return( InvMulConst );
    }
 
    int getLenBits() const
    {
        return( LenBits );
    }
 
    int getLen() const
    {
        return( Len );
    }
 
    void forward( double* const p ) const
    {
    #if R8B_FLOATFFT
 
        float* const op = (float*) p;
        int i;
 
        for( i = 0; i < Len; i++ )
        {
            op[ i ] = (float) p[ i ];
        }
 
    #endif // R8B_FLOATFFT
 
    #if R8B_IPP
 
        ippsFFTFwd_RToPerm_64f( p, p, SPtr, WorkBuffer );
 
    #elif R8B_PFFFT
 
        pffft_transform_ordered( setup, op, op, work, PFFFT_FORWARD );
 
    #elif R8B_PFFFT_DOUBLE
 
        pffftd_transform_ordered( setup, p, p, work, PFFFT_FORWARD );
 
    #else // R8B_PFFFT_DOUBLE
 
        ooura_fft :: rdft( Len, 1, p, wi, wd );
 
    #endif // R8B_IPP
    }
 
    void inverse( double* const p ) const
    {
    #if R8B_IPP
 
        ippsFFTInv_PermToR_64f( p, p, SPtr, WorkBuffer );
 
    #elif R8B_PFFFT
 
        pffft_transform_ordered( setup, (float*) p, (float*) p, work,
            PFFFT_BACKWARD );
 
    #elif R8B_PFFFT_DOUBLE
 
        pffftd_transform_ordered( setup, p, p, work, PFFFT_BACKWARD );
 
    #else // R8B_PFFFT_DOUBLE
 
        ooura_fft :: rdft( Len, -1, p, wi, wd );
 
    #endif // R8B_IPP
 
    #if R8B_FLOATFFT
 
        const float* const ip = (const float*) p;
        int i;
 
        for( i = Len - 1; i >= 0; i-- )
        {
            p[ i ] = ip[ i ];
        }
 
    #endif // R8B_FLOATFFT
    }
 
    void multiplyBlocks( const double* const aip1, const double* const aip2,
        double* const aop ) const
    {
    #if R8B_FLOATFFT
 
        const float* const ip1 = (const float*) aip1;
        const float* const ip2 = (const float*) aip2;
        float* const op = (float*) aop;
 
    #else // R8B_FLOATFFT
 
        const double* const ip1 = aip1;
        const double* const ip2 = aip2;
        double* const op = aop;
 
    #endif // R8B_FLOATFFT
 
    #if R8B_IPP
 
        ippsMulPerm_64f( (Ipp64f*) ip1, (Ipp64f*) ip2, (Ipp64f*) op, Len );
 
    #else // R8B_IPP
 
        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_IPP
    }
 
    void multiplyBlocks( const double* const aip, double* const aop ) const
    {
    #if R8B_FLOATFFT
 
        const float* const ip = (const float*) aip;
        float* const op = (float*) aop;
        float t;
 
    #else // R8B_FLOATFFT
 
        const double* const ip = aip;
        double* const op = aop;
 
        #if !R8B_IPP
        double t;
        #endif // !R8B_IPP
 
    #endif // R8B_FLOATFFT
 
    #if R8B_IPP
 
        ippsMulPerm_64f( (Ipp64f*) op, (Ipp64f*) ip, (Ipp64f*) op, Len );
 
    #else // R8B_IPP
 
        op[ 0 ] *= ip[ 0 ];
        op[ 1 ] *= ip[ 1 ];
 
        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_IPP
    }
 
    void multiplyBlocksZP( const double* const aip, double* const aop ) const
    {
    #if R8B_FLOATFFT
 
        const float* const ip = (const float*) aip;
        float* const op = (float*) aop;
 
    #else // R8B_FLOATFFT
 
        const double* ip = aip;
        double* op = aop;
 
    #endif // R8B_FLOATFFT
 
    // SIMD implementations assume that pointers are address-aligned.
 
    #if !R8B_FLOATFFT && 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 ov3 = _mm_load_pd( op + 4 );
            const __m128d iv4 = _mm_load_pd( ip + 6 );
            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_FLOATFFT && 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( double* const ap ) const
    {
    #if R8B_FLOATFFT
 
        float* const p = (float*) ap;
 
    #else // R8B_FLOATFFT
 
        double* const p = ap;
 
    #endif // R8B_FLOATFFT
 
        int i = 2;
 
        while( i < Len )
        {
            p[ i + 1 ] = p[ i ];
            i += 2;
        }
    }
 
private:
    int LenBits; 
    int Len; 
    double InvMulConst; 
    CDSPRealFFT* Next; 
 
    #if R8B_IPP
        IppsFFTSpec_R_64f* SPtr; 
        CFixedBuffer< unsigned char > SpecBuffer; 
        CFixedBuffer< unsigned char > WorkBuffer; 
    #elif R8B_PFFFT
        PFFFT_Setup* setup; 
        CFixedBuffer< float > work; 
    #elif R8B_PFFFT_DOUBLE
        PFFFTD_Setup* setup; 
        CFixedBuffer< double > work; 
    #else // R8B_PFFFT_DOUBLE
        CFixedBuffer< int > wi; 
        CFixedBuffer< double > wd; 
    #endif // R8B_IPP
 
    class CObjKeeper
    {
        R8BNOCTOR( CObjKeeper );
 
    public:
        CObjKeeper()
            : Object( NULL )
        {
        }
 
        ~CObjKeeper()
        {
            delete Object;
        }
 
        CObjKeeper& operator = ( CDSPRealFFT* const aObject )
        {
            Object = aObject;
            return( *this );
        }
 
        operator CDSPRealFFT* () const
        {
            return( Object );
        }
 
    private:
        CDSPRealFFT* Object; 
    };
 
    CDSPRealFFT()
    {
    }
 
    CDSPRealFFT( const int aLenBits )
        : LenBits( aLenBits )
        , Len( 1 << aLenBits )
    #if R8B_IPP
        , InvMulConst( 1.0 / Len )
    #elif R8B_PFFFT
        , InvMulConst( 1.0 / Len )
    #elif R8B_PFFFT_DOUBLE
        , InvMulConst( 1.0 / Len )
    #else // R8B_PFFFT_DOUBLE
        , InvMulConst( 2.0 / Len )
    #endif // R8B_IPP
    {
    #if R8B_IPP
 
        int SpecSize;
        int SpecBufferSize;
        int BufferSize;
 
        ippsFFTGetSize_R_64f( LenBits, IPP_FFT_NODIV_BY_ANY,
            ippAlgHintFast, &SpecSize, &SpecBufferSize, &BufferSize );
 
        CFixedBuffer< unsigned char > InitBuffer( SpecBufferSize );
        SpecBuffer.alloc( SpecSize );
        WorkBuffer.alloc( BufferSize );
 
        ippsFFTInit_R_64f( &SPtr, LenBits, IPP_FFT_NODIV_BY_ANY,
            ippAlgHintFast, SpecBuffer, InitBuffer );
 
    #elif R8B_PFFFT
 
        setup = pffft_new_setup( Len, PFFFT_REAL );
        work.alloc( Len );
 
    #elif R8B_PFFFT_DOUBLE
 
        setup = pffftd_new_setup( Len, PFFFT_REAL );
        work.alloc( Len );
 
    #else // R8B_PFFFT_DOUBLE
 
        wi.alloc( (int) ceil( 2.0 + sqrt( (double) ( Len >> 1 ))));
        wi[ 0 ] = 0;
        wd.alloc( Len >> 1 );
 
    #endif // R8B_IPP
    }
 
    ~CDSPRealFFT()
    {
        #if R8B_PFFFT
            pffft_destroy_setup( setup );
        #elif R8B_PFFFT_DOUBLE
            pffftd_destroy_setup( setup );
        #endif // R8B_PFFFT_DOUBLE
 
        delete Next;
    }
};
 
class CDSPRealFFTKeeper : public R8B_BASECLASS
{
    R8BNOCTOR( CDSPRealFFTKeeper );
 
public:
    CDSPRealFFTKeeper()
        : Object( NULL )
    {
    }
 
    CDSPRealFFTKeeper( const int LenBits )
    {
        Object = acquire( LenBits );
    }
 
    ~CDSPRealFFTKeeper()
    {
        if( Object != NULL )
        {
            release( Object );
        }
    }
 
    const CDSPRealFFT* operator -> () const
    {
        R8BASSERT( Object != NULL );
 
        return( Object );
    }
 
    void init( const int LenBits )
    {
        if( Object != NULL )
        {
            if( Object -> LenBits == LenBits )
            {
                return;
            }
 
            release( Object );
        }
 
        Object = acquire( LenBits );
    }
 
    void reset()
    {
        if( Object != NULL )
        {
            release( Object );
            Object = NULL;
        }
    }
 
private:
    CDSPRealFFT* Object; 
 
    static CSyncObject StateSync; 
    static CDSPRealFFT :: CObjKeeper FFTObjects[]; 
 
    CDSPRealFFT* acquire( const int LenBits )
    {
        R8BASSERT( LenBits > 0 && LenBits <= 30 );
 
        R8BSYNC( StateSync );
 
        if( FFTObjects[ LenBits ] == NULL )
        {
            return( new CDSPRealFFT( LenBits ));
        }
 
        CDSPRealFFT* ffto = FFTObjects[ LenBits ];
        FFTObjects[ LenBits ] = ffto -> Next;
 
        return( ffto );
    }
 
    void release( CDSPRealFFT* const ffto )
    {
        R8BSYNC( StateSync );
 
        ffto -> Next = FFTObjects[ ffto -> LenBits ];
        FFTObjects[ ffto -> LenBits ] = ffto;
    }
};
 
inline void calcMinPhaseTransform( double* const Kernel, const int KernelLen,
    const int LenMult = 2, const bool DoFinalMul = true,
    double* const DCGroupDelay = 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< double > ip2( Len2 + 1 );
 
    memcpy( &ip[ 0 ], Kernel, KernelLen * sizeof( ip[ 0 ]));
    memset( &ip[ KernelLen ], 0, ( Len - KernelLen ) * sizeof( ip[ 0 ]));
 
    CDSPRealFFTKeeper ffto( LenBits );
    ffto -> forward( ip );
 
    // Create the "log |c|" spectrum while saving the original power spectrum
    // in the "ip2" buffer.
 
    #if R8B_FLOATFFT
        float* const aip = (float*) &ip[ 0 ];
        float* const aip2 = (float*) &ip2[ 0 ];
        const float nzbias = 1e-35;
    #else // R8B_FLOATFFT
        double* const aip = &ip[ 0 ];
        double* const aip2 = &ip2[ 0 ];
        const double nzbias = 1e-300;
    #endif // R8B_FLOATFFT
 
    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( 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.
 
    ffto -> forward( ip );
 
    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( ip );
 
    if( DoFinalMul )
    {
        for( i = 0; i < KernelLen; i++ )
        {
            Kernel[ i ] = ip[ i ] * m1;
        }
    }
    else
    {
        memcpy( &Kernel[ 0 ], &ip[ 0 ], KernelLen * sizeof( Kernel[ 0 ]));
    }
 
    if( DCGroupDelay != NULL )
    {
        *DCGroupDelay = calcFIRFilterGroupDelay( Kernel, KernelLen, 0.0 );
    }
}
 
} // namespace r8b
 
#endif // VOX_CDSPREALFFT_INCLUDED