viterbi.cpp

来自「signal-processing.rar信号处理demo原码」· C++ 代码 · 共 562 行 · 第 1/2 页

        viterbiInf->ceInput = *ce32Input;
        break;

    case ViterbiState64:
        viterbiInf->vState = 64;
        viterbiInf->numPossStates = 16;
        viterbiInf->ceStates = *ceStates64;
        viterbiInf->ceInput = *ce64Input;
        break;

    default:
        break;
    }

    viterbiInf->ceDlState = 0;

    return;
} // ViterbiInit()


bool TrellisEncoder( const Real4DSymbl *mapOut, ViterbiTestInfo *vInf )
{
    bool    U0;                                 /* trellis encoder output bit           */
    Ipp8u   ceInput;                            /* convolutional encoder input Y4Y3Y2Y1 */
    Ipp8u   ssLabel_0, ssLabel_1;               /* 2 subset labels s(2m), s(2m+1)       */
    Ipp8u   *ceState = &vInf->ceDlState;        /* convolutional encoder delay state    */

    U0 = *ceState & 0x1;

/* --------------------- Symbol-to-bit converter --------------------------- */
    ssLabel_0 = subsetLabel[(35 - mapOut->m_firstPoint.im)/2][(35 + mapOut->m_firstPoint.re)/2];
    ssLabel_1 = subsetLabel[(35 - mapOut->m_secondPoint.im)/2][(35 + mapOut->m_secondPoint.re)/2];

/* --------------------- Convolutional encoder ----------------------------- */
    ceInput = vInf->ceInput[8*ssLabel_0 + ssLabel_1];
    *ceState = vInf->ceStates[ *ceState*vInf->numPossStates + ceInput];
    U0 = *ceState & 0x1;

    return U0;
} // TrellisEncoder()


void GetTrellisSequence( Real4DSymbl *trSeq, int len, ViterbiTestInfo *vInf )
{
    int         i;
    Ipp8u       I1;     /* bit I1(i,j) from parser                          */
    Ipp8u       U0;     /* trellis encoder output bit                       */
    Ipp8u       Z;      /* differential encoder output integer              */
    Real4DSymbl v;      /* 2 2D points from the quarter superconstellation  */
    Real4DSymbl u;      /* 2 2D output siganal points                       */

/* ------------------ Quarter superconstellation indices ------------------- */
    const Ipp16s scIdx[18] = { -35, -31, -27, -23, -19, -15, -11, -7, -3,
                                 1,   5,   9,  13,  17,  21,  25, 29, 33 };
    U0 = 0;

    srand((unsigned int) time(NULL) );

    for( i = 0; i < len; i++ )
    {
/* --------------------- Initialize input data ----------------------------- */
        v.m_firstPoint.re = scIdx[rand()/1821];
        v.m_firstPoint.im = scIdx[rand()/1821];

        v.m_secondPoint.re = scIdx[rand()/1821];
        v.m_secondPoint.im = scIdx[rand()/1821];

        I1 = (Ipp8u)(rand()/16384);
      
        Z = (Ipp8u)(rand()/10923);
       
/* ---------------- Rotate it according to Z, I1, U0 ----------------------- */
        switch(Z)   /* rotate first point  */
        {
        case 0:      /* 0 degrees           */
            u.m_firstPoint.re = v.m_firstPoint.re;
            u.m_firstPoint.im = v.m_firstPoint.im;
            break;

        case 1:     /* 90 degrees         */
            u.m_firstPoint.re = v.m_firstPoint.im;
            u.m_firstPoint.im = (Ipp16u)(-v.m_firstPoint.re);
            break;

        case 2:     /* 180 degrees        */
            u.m_firstPoint.re = (Ipp16u)(-v.m_firstPoint.re);
            u.m_firstPoint.im = (Ipp16u)(-v.m_firstPoint.im);
            break;

        case 3:    /* 270 degrees        */
            u.m_firstPoint.re = (Ipp16u)(-v.m_firstPoint.im);
            u.m_firstPoint.im = v.m_firstPoint.re;
            break;

        default:
            u.m_firstPoint.re = 0;
            u.m_firstPoint.im = 0;
        }

        switch( (Z + 2*I1 + U0)%4 ) /* rotate second point  */
        {
        case 0:      /* 0 degrees         */
            u.m_secondPoint.re = v.m_secondPoint.re;
            u.m_secondPoint.im = v.m_secondPoint.im;
            break;

        case 1:     /* 90 degrees        */
            u.m_secondPoint.re = v.m_secondPoint.im;
            u.m_secondPoint.im = (Ipp16u)(-v.m_secondPoint.re);
            break;

        case 2:     /* 180 degrees       */
            u.m_secondPoint.re = (Ipp16u)(-v.m_secondPoint.re);
            u.m_secondPoint.im = (Ipp16u)(-v.m_secondPoint.im);
            break;

        case 3:     /* 270 degrees       */
            u.m_secondPoint.re = (Ipp16u)(-v.m_secondPoint.im);
            u.m_secondPoint.im = v.m_secondPoint.re;
            break;

        default:
            u.m_secondPoint.re = 0;
            u.m_secondPoint.im = 0;
       }

/* -------------------------------- Output --------------------------------- */
        trSeq[i].m_firstPoint.re = u.m_firstPoint.re;
        trSeq[i].m_firstPoint.im = u.m_firstPoint.im;
        trSeq[i].m_secondPoint.re = u.m_secondPoint.re;
        trSeq[i].m_secondPoint.im = u.m_secondPoint.im;

/* -------------------------------- Get U0 --------------------------------- */
        U0 = TrellisEncoder( &u, vInf );
    }

    return;
} // GetTrellisSequence()


int AddChannelNoise( Real4DSymbl *tOutput, Real4DSymbl *vInput,int len, double SNR )
{
    int             i, j;
    int             sLen;           /* signal lenght            */
    float           sMagn;          /* output signal power      */
    float           cFreq;          /* carrier frequency        */
    float           pPhase;         /* tone phase               */
    float           *samples;       /* output signal            */
    float           *Cos, *Sin;
    Ipp16s          *gns;           /* white gaussian nois      */
    Ipp16s          stDev;          /* noise standard deviation */
    unsigned int    gnsSeed;
    IppStatus       st;
    sLen = 6*len;

/* ----------------------- Generate output samples ------------------------- */
    Cos = ippsMalloc_32f( sLen );
    
    Sin = ippsMalloc_32f( sLen );
    
    cFreq = 0.1904345f;    /* maximum speed */

    pPhase = 0.0f;
    st = ippsTone_Direct_32f( Cos, sLen, 1.0f, cFreq, &pPhase, ippAlgHintAccurate );
    if(ippStsNoErr != st)
    {
       printf(" Error in ippsToneDirect: %s\n", ippGetStatusString(st));
       return 1;
   }
    pPhase = (float)IPP_PI2;
    st = ippsTone_Direct_32f( Sin, sLen, 1.0f, cFreq, &pPhase, ippAlgHintAccurate );
    if(ippStsNoErr != st)
    {
       printf(" Error in ippsTone_Direct_32f: %s\n", ippGetStatusString(st));
       return 1;
    }
    samples = ippsMalloc_32f( sLen );

    for( i = 0; i < len; i++ )
    {
        tOutput[i].m_firstPoint.re = (Ipp16s)(128*tOutput[i].m_firstPoint.re);
        tOutput[i].m_firstPoint.im = (Ipp16s)(128*tOutput[i].m_firstPoint.im);

        for( j = 6*i; j < 6*i+3; j++ )
        {
            samples[j] = tOutput[i].m_firstPoint.re * Cos[j] - tOutput[i].m_firstPoint.im * Sin[j];
        }

        tOutput[i].m_secondPoint.re = (Ipp16s)(128*tOutput[i].m_secondPoint.re);
        tOutput[i].m_secondPoint.im = (Ipp16s)(128*tOutput[i].m_secondPoint.im);

        for( j = 6*i+3; j < 6*i+6; j++ )
        {
            samples[j] = tOutput[i].m_secondPoint.re * Cos[j] - tOutput[i].m_secondPoint.im * Sin[j];
        }
    }

    ippsFree( Cos );
    ippsFree( Sin );

/* ---------------------- Compute output signal power ---------------------- */
    st = ippsNorm_L2_32f( samples, sLen, &sMagn );
    if(ippStsNoErr != st)
    {
       printf(" Error in ippsNorm_L2_32f: %s\n", ippGetStatusString(st));
       return 1;
    }
    ippsFree( samples );

/* ------------------- Compute noise standard deviation -------------------- */
    stDev = (Ipp16s)( sMagn/sqrt(sLen*pow(10.0, SNR/10.0)) );

/* ------------------- Generate white gaussian noise ----------------------- */
    gns = ippsMalloc_16s( sLen );
   
    st = ippsRandGauss_Direct_16s( gns, sLen, 0, stDev, &gnsSeed );
    if(ippStsNoErr != st)
    {
       printf(" Error in ippsRandGauss_Direct_16s: %s\n", ippGetStatusString(st));
       return 1;
    }
/* ---------------- Add noise to Viterbi decoder input --------------------- */
    for( i = 0; i < len; i++ )
    {
        vInput[i].m_firstPoint.re  = (Ipp16s)(tOutput[i].m_firstPoint.re + gns[4*i]);
        vInput[i].m_firstPoint.im  = (Ipp16s)(tOutput[i].m_firstPoint.im + gns[4*i+1]);
        vInput[i].m_secondPoint.re = (Ipp16s)(tOutput[i].m_secondPoint.re + gns[4*i+2]);
        vInput[i].m_secondPoint.im = (Ipp16s)(tOutput[i].m_secondPoint.im + gns[4*i+3]);
    }

    ippsFree( gns );

    return 0;
} // AddChannelNoise()


int CheckViterbiOut( const Real4DSymbl     *tOutput,
                     const Real4DSymbl     *vInput,
                     int                   len,
                     const ViterbiTestInfo *vInf,
                     long                  *failCount)
{
    int         i;
    Ipp16sc     inputPoint[2];          /* 4D Viterbi decoder input point   */
    Ipp16sc     outputPoint[2];         /* 4D Viterbi decoder output point  */
    ViterbiInfo viterbiParam;
   
    viterbiParam.m_State = vInf->vState;

/* ------------------------ Initialise Viterbi decoder --------------------- */
    InitViterbiDecoder( &viterbiParam);

/* --------------------------- Run Viterbi decoder ------------------------- */
    for( i = 0; i < len; i++ )
    {
        inputPoint[0].re = vInput[i].m_firstPoint.re;
        inputPoint[0].im = vInput[i].m_firstPoint.im;
        inputPoint[1].re = vInput[i].m_secondPoint.re;
        inputPoint[1].im = vInput[i].m_secondPoint.im;
        
        if(GetViterbiData( inputPoint, outputPoint, &viterbiParam))
        {
            return 1;
        }

        if( i >= viterbiParam.m_State )
        {
/* -------------------------- Check Viterbi output ------------------------- */
            if( outputPoint[0].re != tOutput[i-viterbiParam.m_State].m_firstPoint.re  ||
                outputPoint[0].im != tOutput[i-viterbiParam.m_State].m_firstPoint.im  ||
                outputPoint[1].re != tOutput[i-viterbiParam.m_State].m_secondPoint.re ||
                outputPoint[1].im != tOutput[i-viterbiParam.m_State].m_secondPoint.im   )
            {   
                (*failCount)++;
            }
        }
    }

    return 0;
} // CheckViterbiOut()