transform.h

/******** Digital signal processing library v0.01 ********
                
                Header file Transform.h
        Contains interfaces for abstract class ITransform, 
        Template classes: WFastHT, WFastFT, WFastWT
                                     WIFastHT
                                    
                Template class WFastHT - implements algorighm of 
        Fast Haara's Transformation (by Anatoly Beletsky). Template has to type parameters:
                T_in            - type of input values (input signal)
                T_out   - type of output values (spectrum).     
        Input signal has any type(int,double, complex<T> ).
        Spectrum is double or complex.
        
                Template class WIFastHT - implements algorithm of
        Inverse Fast Haara's Transformation.
        
                Template class WFastFT - implements algorithm of
        Fast Fourier Transformation (Tim Kientzle "A Programmers Guide To Sound"
        Copyright (c) 1998 Tim Kientzle)
        
                Template class WFastWT - implements algorihtm of
        Fast Walsh's Transformation in different basises.
                HADOMARD - Walshs-Hadomars basis
                PELEY      - Peleys basis
                WALSH      - "classical" basis
                COLEY      - Coleys basis (discovered by Anatoly Belestky in 1999)

Alexander Beletsky
e-mail:  albel@pochtamt.ru

À.À. Beletsky 2002(c)
*********************************************************/

#ifndef Transform_H
#define Transform_H

#include <cmath>
#include <complex>
using namespace std;


#include "transerr.h"

const double PI = 3.14159265358979323846264338;
//Basises of Fast Walshs Transformation
enum WBasis {
        HADOMARD,       
        PELEY,          
        WALSH,          
        COOLEY
};

//--------------------------------------------------------------

template < class T_in, class T_out >
class ITransform {
public:
        virtual void setData( T_in *dat ) = 0;
        virtual void doTransform() = 0;
};

//--------------------------------------------------------------

// Template class WFastHT
template< class T_in, class T_out >
class WFastHT : public ITransform< T_in, T_out > {
public:
        WFastHT( uint N_par = 8 );
        void setData( T_in *dat );
        void doTransform();
        T_out *getSpectrum();
        virtual ~WFastHT();
private:
        void regroupe();
        void multiply();
        uint binInv( uint index );
        T_out *spectrum;
        T_out *aux;
        uint N, m;
        bool noData;
};

template< class T_in, class T_out >
WFastHT< T_in, T_out >::WFastHT( uint N_par ) : N( N_par ), noData( true )
{
        //is N binary rational?
        if( ((~N + 1)&N) != N ) {
                throw TransError::badNValue();
        }
        //Log(n) - computation
        int count = 0;
        for(int check = N_par; (check&1) == 0 ; count++, check >>= 1 );
        m = count;

        spectrum = new T_out[N];
        aux = new T_out[N];
}

template< class T_in, class T_out >
WFastHT< T_in, T_out >::~WFastHT()
{
        delete [] spectrum;
        delete [] aux;
}

template< class T_in, class T_out >
void WFastHT< T_in, T_out >::setData( T_in *dat )
{
        if( dat != 0 ) {
                for(int i = 0; i < N; i++ ) aux[i] = static_cast< T_out >( dat[i] );
                noData = false;
        }
        else throw TransError::noData();
}

template< class T_in, class T_out >
void WFastHT< T_in, T_out >::doTransform()
{
        if( noData )
                throw TransError::noData();
        else {
        //Spectrum computation
        for( int stage = 0, make = N/2, turn = 1; stage < m; stage++, make /= 2, turn *= 2 ) {
                for( int step = 0; step < make; step++ ) {
                        //"butterfly" operation
                        T_out A = aux[2*turn*step] - aux[2*turn*step + turn];           
                        aux[2*turn*step] += aux[2*turn*step + turn];
                        aux[2*turn*step + turn] = A;
                }
        }
        //Regrouping of harmonics
        regroupe();
        //Multiplication with nessesary coeffisients
        multiply();
        }
}

template< class T_in, class T_out >
T_out *WFastHT< T_in, T_out>::getSpectrum()
{
        return spectrum;
}

template< class T_in, class T_out >
void WFastHT< T_in, T_out >::regroupe() 
{
        //Regroping of harmonics

        for( uint i = 0, k = 0; i < N; i++ ) {
                if( i < N/2 ) {
                        spectrum[i] = aux[ binInv( i ) ];
                }
                else {
                        spectrum[i] = aux[2*k + 1];
                        k++;
                }
        }
}

template< class T_in, class T_out >
void WFastHT< T_in, T_out >::multiply()
{
        spectrum[0] /= N;
        spectrum[1] /= N;
        
        double S = sqrt(2.0), coeff = 1;
        for(int i = 1, l_1 = 2, l_2 = l_1*2; i < m; i++, l_1 *= 2, l_2 = l_1*2 ) {
                coeff *= S;
                for(int j = l_1; j < l_2; j++) {
                        spectrum[j] *= coeff;
                        spectrum[j] /= N;
                }
        }
}

//Binary inversion operation
template< class T_in, class T_out >
uint WFastHT< T_in, T_out >::binInv( uint index )
{
        uint LMASK = 1, RMASK = N >> 1, inverse = index;
        bool nextIter = true;
        
        while( nextIter ) {
                if( LMASK == RMASK ) nextIter = false;
                if( ( (index & LMASK) == 0 ) && ( (index & RMASK) == 0) ||
                        ( (index & LMASK) != 0 ) && ( (index & RMASK) != 0) ) {
                //return inverse;
                }
                else {
                        if( (index & LMASK) == LMASK ) {
                                inverse |= RMASK;
                                inverse &= ~LMASK;
                        }
                        else {
                                if( (index & RMASK) == RMASK ) {
                                        inverse |= LMASK;
                                        inverse &= ~RMASK;
                                }
                        }
                }
                LMASK <<= 1;
                RMASK >>= 1;
        }
        return inverse;
}

//--------------------------------------------------------------

//Template class WIFastFT
template< class T_in >
class WFastFT : public ITransform< T_in, complex< double > > {
public:
        WFastFT( uint N_par = 8 );
        void setData( T_in *dat );
        void doTransform();
        complex< double > *getSpectrum();
        virtual ~WFastFT();
private:
        void regroupe();
        uint binInv( uint index ); 
        T_in *inputData;
        complex< double > *spectrum;
        uint N;
        bool noData;
};

template< class T_in >
WFastFT< T_in >::WFastFT( uint N_par ) : N( N_par ), noData( true )
{
        if( ((~N + 1)&N) != N ) {
                throw TransError::badNValue();
        }

        spectrum = new complex< double >[N];
}

template< class T_in >
WFastFT< T_in >::~WFastFT()
{
        delete [] spectrum;
}

template< class T_in >
void WFastFT< T_in >::setData( T_in *dat )
{
        if( dat != 0 ) {
                inputData = dat;

                regroupe();
                noData = false;
        }
        else throw TransError::noData();
}

template< class T_in >
void WFastFT< T_in >::doTransform()
{
        if( noData )
                throw TransError::noData();
        else {
                for( uint halfSize = 1; halfSize < N; halfSize *= 2) {
                        complex< double > phaseShiftStep = polar( 1.0 , -PI/halfSize );
                        complex< double > currentPhaseShift = polar( 1.0 );
                        for( uint fftStep = 0; fftStep < halfSize; fftStep++ ) {
                                for( uint i = fftStep; i < N; i += 2*halfSize ) {
                                        complex< double > A = currentPhaseShift * spectrum[i + halfSize];
                                        spectrum[i + halfSize] = spectrum[i] - A;
                                        spectrum[i] += A;
                                }
                                currentPhaseShift *= phaseShiftStep;
                        }
                }
        }
}

template< class T_in >
complex< double > *WFastFT< T_in >::getSpectrum()
{
        return spectrum;
}

template< class T_in >
void WFastFT< T_in >::regroupe()
{
        for( uint i = 0; i < N; i++ ) {
                spectrum[i] = static_cast< complex< double > >( inputData[ binInv( i ) ] );
        }
}

template< class T_in >
uint WFastFT< T_in >::binInv( uint index  )
{
        uint LMASK = 1, RMASK = N >> 1, inverse = index;
        bool nextIter = true;
        
        while( nextIter ) {
                if( LMASK == RMASK ) nextIter = false;
                if( ( (index & LMASK) == 0 ) && ( (index & RMASK) == 0) ||
                        ( (index & LMASK) != 0 ) && ( (index & RMASK) != 0) ) {
                //return inverse;
                }
                else {
                        if( (index & LMASK) == LMASK ) {
                                inverse |= RMASK;
                                inverse &= ~LMASK;
                        }
                        else {
                                if( (index & RMASK) == RMASK ) {
                                        inverse |= LMASK;
                                        inverse &= ~RMASK;
                                }
                        }
                }
                LMASK <<= 1;
                RMASK >>= 1;
        }
        return inverse;
}

//--------------------------------------------------------------

//Template class WFastWT
template < class T_in, class T_out >
class WFastWT : ITransform< T_in, T_out > {
public:
        WFastWT( uint N_par = 8, WBasis bas = HADOMARD );
        void setData( T_in *dat );
        void doTransform();
        T_out *getSpectrum();
        virtual ~WFastWT();
private:
        void regroupe();
        uint binInv( uint index );
        uint binInvGCL( uint index );
        uint binInvGCR( uint index );
        T_in *inputData;
        T_out *spectrum;
        uint N;
        WBasis basis;
        bool noData;
};

template< class T_in, class T_out >
WFastWT< T_in, T_out>::WFastWT( uint N_par, WBasis bas ) : N( N_par ), basis( bas ),
                                                                                                                        noData( true )
{
        if( ((~N + 1)&N) != N ) {
                throw TransError::badNValue();
        }

        spectrum = new T_out[N];
}

template< class T_in, class T_out >
WFastWT< T_in, T_out >::~WFastWT()
{
        delete [] spectrum;
}

template< class T_in, class T_out >
void WFastWT< T_in, T_out >::setData( T_in *dat )
{
        if( dat != 0 ) {
                inputData = dat;

                regroupe();
                noData = false;
        }
        else throw TransError::noData();
}

template< class T_in, class T_out >
void WFastWT< T_in, T_out >::doTransform()
{
        if( noData ) {
                throw TransError::noData();
        }
        else {
                for( int halfSize = 1; halfSize < N; halfSize *= 2) {
                        for( int fwtStep = 0; fwtStep < halfSize; fwtStep++ ) {
                                for( int i = fwtStep; i < N; i += 2*halfSize ) {
                                        T_out A = spectrum[i + halfSize];
                                        spectrum[i + halfSize] = spectrum[i] - A;
                                        spectrum[i] += A;
                                }
                        }
                }
        }
}

template< class T_in, class T_out >
T_out *WFastWT< T_in, T_out >::getSpectrum()
{
        return spectrum;
}

template< class T_in, class T_out >
void WFastWT< T_in, T_out >::regroupe()
{
        int i;
        switch( basis ) {
        case HADOMARD:
                for( i = 0; i < N; i++ ) {
                        spectrum[i] = static_cast< T_out >( inputData[i] );
                }
                break;
        case PELEY:
                for( i = 0; i < N; i++ ) {
                        spectrum[i] = static_cast< T_out >( inputData[ binInv( i ) ] );
                }
                break;
        case WALSH:
                for( i = 0; i < N; i++ ) {
                        int newI = binInv( i );
                        spectrum[i] = static_cast< T_out >( inputData[ binInvGCL( newI ) ] );
                }
                break;
        case COOLEY:
                for( i = 0; i < N; i++ ) {
                        int newI = binInv( i );
                        spectrum[i] = static_cast< T_out >( inputData[ binInvGCR( newI ) ] );
                }
                break;
        default:
                throw TransError::noSuchBasis();
        }
}

template< class T_in, class T_out >
uint WFastWT< T_in, T_out >::binInv( uint index  )
{
        //base >>= 1;
        uint LMASK = 1, RMASK = N >> 1, inverse = index;
        bool nextIter = true;
        
        while( nextIter ) {
                if( LMASK == RMASK ) nextIter = false;
                if( ( (index & LMASK) == 0 ) && ( (index & RMASK) == 0) ||
                        ( (index & LMASK) != 0 ) && ( (index & RMASK) != 0) ) {
                //return inverse;
                }
                else {
                        if( (index & LMASK) == LMASK ) {
                                inverse |= RMASK;
                                inverse &= ~LMASK;
                        }
                        else {
                                if( (index & RMASK) == RMASK ) {
                                        inverse |= LMASK;
                                        inverse &= ~RMASK;
                                }
                        }
                }
                LMASK <<= 1;
                RMASK >>= 1;
        }
        return inverse;
}

//Binary operation of Greys lefthanded transposition (classical coding)
template< class T_in, class T_out >
uint WFastWT< T_in, T_out >::binInvGCL( uint index  )
{
        unsigned int LMASK = N >> 1, RMASK = LMASK >> 1, greyCL = 0;

        if( ( index & LMASK ) != 0 )
                 greyCL |= LMASK;

        while( RMASK ) {
                if( ( ( greyCL & LMASK ) != 0 ) && ( ( index & RMASK ) == 0 ) ||
                        ( ( greyCL & LMASK ) == 0 ) && ( ( index & RMASK ) != 0 ) ) {
                        greyCL |= RMASK;
                }
                LMASK >>= 1;
                RMASK >>= 1;
        }
        return greyCL;
}

//Binary operation of Greys righthanded transposition (discovered by Anatoly Beletsky)
template< class T_in, class T_out >
uint WFastWT< T_in, T_out >::binInvGCR( uint index  )
{
        unsigned int base = N >> 1;
        unsigned int LMASK = 1, RMASK = LMASK << 1, greyCR = 0;

        if( ( index & LMASK ) != 0 )
                greyCR |= LMASK;

        while( LMASK != base ) {
                if( ( ( greyCR & LMASK ) != 0 ) && ( ( index & RMASK ) == 0 ) ||
                    ( ( greyCR & LMASK ) == 0 ) && ( ( index & RMASK ) != 0 ) ) {
                        greyCR |= RMASK;
                }
                LMASK <<= 1;
                RMASK <<= 1;
        }
        return greyCR;
}

//-------------------------------------------------------------------

template < class T_in, class T_out >
class WIFastHT : public ITransform< T_in, T_out > {
public:
        WIFastHT( uint N_par = 8 );
        void doTransform();
        void setData( T_in *dat );
        T_out *getSpectrum();
        virtual void regroupe();
        virtual ~WIFastHT();
private:
        uint N, m;
        void multiply();
        T_out *spectrum;
        T_out *aux1;
        T_out *aux2;
        bool noData;
};

template< class T_in, class T_out >
WIFastHT< T_in, T_out>::WIFastHT( uint N_par ) : N( N_par ), noData( true )
{

        if( ((~N + 1)&N) != N ) {
                throw TransError::badNValue();
        }

        int count = 0;
        for(int check = N_par; (check&1) == 0 ; count++, check >>= 1 );
        m = count;

        spectrum = new T_out[ N ];
        aux1 = new T_out[ N ];
        aux2 = new T_out[ N ];
}

template< class T_in, class T_out >
WIFastHT< T_in, T_out >::~WIFastHT()
{
        delete [] spectrum;
        delete [] aux1;
        delete [] aux2;
}

template< class T_in, class T_out >
void WIFastHT< T_in, T_out >::doTransform()
{
        if( noData )
                throw TransError::noData();
        else {
                spectrum[ 0 ] = aux1[ 0 ] + aux1[ 1 ];
                spectrum[ 1 ] = aux1[ 0 ] - aux1[ 1 ];
                for( int halfSize = 2; halfSize < N; halfSize *= 2) {
                        for( int i = 0; i < halfSize; i++ ) aux2[ i ] = spectrum[ i ]; 
                        uint count = halfSize;
                        for( int ifhtStep = 0; ifhtStep < halfSize; ifhtStep++ ) {
                                spectrum[ 2*ifhtStep ] = aux2[ ifhtStep ] + aux1[ count ];
                                spectrum[ 2*ifhtStep + 1 ] = aux2[ ifhtStep ] - aux1[ count ];
                                count++;
                        }
                }
        }
}

template< class T_in, class T_out >
void WIFastHT< T_in, T_out >::setData( T_in *dat )
{
        for( int i = 0; i < N; i++ ) {
                aux1[i] = static_cast< T_out >( dat[i] );
        }
        multiply();
        noData = false;
}

template< class T_in, class T_out >
void WIFastHT< T_in, T_out >::multiply()
{
        double S = sqrt(2.0), coeff = 1;
        for(int i = 1, l_1 = 2, l_2 = l_1*2; i < m; i++, l_1 *= 2, l_2 = l_1*2 ) {
                coeff *= S;
                for(int j = l_1; j < l_2; j++) {
                        aux1[j] *= coeff;
                }
        }
}

template< class T_in, class T_out >
T_out *WIFastHT< T_in, T_out >::getSpectrum()
{
        return spectrum;
}

template< class T_in, class T_out >
void WIFastHT< T_in, T_out >::regroupe()
{

}

#endif

---