CDSPRealFFT.h

Go to the documentation of this file.

//$ nobt
//$ nocpp

#ifndef R8B_CDSPREALFFT_INCLUDED
#define R8B_CDSPREALFFT_INCLUDED

#include "r8bbase.h"

#if !R8B_IPP && !R8B_PFFFT
    #include "fft4g.h"
#endif // !R8B_IPP && !R8B_PFFFT

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 );

    #else // R8B_PFFFT

        ooura_fft :: rdft( Len, 1, p, wi.getPtr(), wd.getPtr() );

    #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 );

    #else // R8B_PFFFT

        ooura_fft :: rdft( Len, -1, p, wi.getPtr(), wd.getPtr() );

    #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;
        double t;

    #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 multiplyBlocksZ( 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* const ip = aip;
        double* const op = aop;

    #endif // R8B_FLOATFFT

    #if R8B_IPP

        ippsMul_64f_I( (const Ipp64f*) ip, (Ipp64f*) op, Len );

    #else // R8B_IPP

        int i;

        for( i = 0; i < Len; i++ )
        {
            op[ i ] *= ip[ i ];
        }

    #endif // R8B_IPP
    }

    void convertToZ( 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; 
    #else // R8B_PFFFT
        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 )
    #else // R8B_PFFFT
        , 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 );

    #else // R8B_PFFFT

        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 );
        #endif // R8B_PFFFT

        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( double ));
    memset( &ip[ KernelLen ], 0, ( Len - KernelLen ) * sizeof( double ));

    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( double ));
    }

    if( DCGroupDelay != NULL )
    {
        double tmp;

        calcFIRFilterResponseAndGroupDelay( Kernel, KernelLen, 0.0,
            tmp, tmp, *DCGroupDelay );
    }
}

} // namespace r8b

#endif // VOX_CDSPREALFFT_INCLUDED