turbo.cpp

基于VC++的turbo码MAP译码算法

            b[x]=NULL;
        }
    }

    if (initialized == false || size != lastsize)
    {
        initialized=true;
        lastsize=size;


        // create the arrays dynamically at runtime, delete at end of routine
        // create gamma[encoder][uk][k][state]
        for (int y=0;y<2;y++)
            for (int x=0;x<2;x++)
            {
                gamma[y][x]=new double*[size];
                gammaE[y][x]=new double*[size];
                for (int z=0;z<size;z++)
                {
                    gamma[y][x][z]=new double[numstates];
                    gammaE[y][x][z]=new double[numstates];
                }
            }

        // create L[encoder #]
        for (int y=0;y<2;y++)
            L[y]=new double[size];

        // create alpha[encoder #][k][state]
        for (int x=0;x<2;x++)
        {
            a[x]=new double*[size];
            b[x]=new double*[size];
            // each Yk has 'numstates' values of gamma
            for (int y=0;y<size;y++)
            {
                a[x][y]=new double[numstates];
                b[x][y]=new double[numstates];
            }
        }
    }

    // initialization of iteration arrays
    for (int x=0;x<numstates;x++)
        a[0][0][x]=b[0][size-1][x]=a[1][0][x]=(x==0) ? 1.0 : 0.0;

    // initialization of extrinsic information array from decoder 2,
    // used in decoder 1
    for (int x=0;x<size;x++)
        L[1][x]=0.0;

    // 4*Eb/No
    // double Lc=4.0/No;
    double Lc=2.0/sigma;

    for (int c=0;c<ITERATIONS;c++)
    {
        // k from 0 to N-1 instead of 1 to N
        for (int k=0;k<size;k++)
        {
            // calculate the gamma(s',s);
            for (int input=0;input<2;input++)
            {
                double uk=(input == 0) ? -1.0 : 1.0;

                for (int s=0;s<numstates;s++)
                {
                    double xk=(output[input][s] == 0) ? -1.0 : 1.0;

                    gammaE[0][input][k][s]=exp(0.5*Lc*parity1[k]*xk);
                    gamma[0][input][k][s] =exp(0.5*uk*(L[1][k]+Lc*mesg[k]))*gammaE[0][input][k][s];
                }
            }
        }

        // calculate the alpha terms
        // from 1 to N-1, 0 is precalculated, N is never used
        for (int k=1;k<size;k++)
        {
            double temp=0;
            // calculate denominator
            for (int state=0;state<numstates;state++)
                temp+=a[0][k-1][fromstate[0][state]]*gamma[0][0][k-1][fromstate[0][state]]+a[0][k-1][fromstate[1][state]]*gamma[0][1][k-1][fromstate[1][state]];
            for (int state=0;state<numstates;state++)
                a[0][k][state]=( a[0][k-1][fromstate[0][state]]*gamma[0][0][k-1][fromstate[0][state]]+a[0][k-1][fromstate[1][state]]*gamma[0][1][k-1][fromstate[1][state]])/temp;
        }

        // calculate the beta terms
        // from N-1 to 1
        for (int k=size-1;k>=1;k--)
        {
            double temp=0;
            // calculate denominator
            for (int state=0;state<numstates;state++)
                temp += a[0][k][fromstate[0][state]]*gamma[0][0][k][fromstate[0][state]]+a[0][k][fromstate[1][state]]*gamma[0][1][k][fromstate[1][state]];
            for (int state=0;state<numstates;state++)
                b[0][k-1][state]=( b[0][k][tostate[0][state]]*gamma[0][0][k][state]+b[0][k][tostate[1][state]]*gamma[0][1][k][state])/temp;
        }

        // compute the extrinsic LLRs
        for (int k=0;k<size;k++)
        {
            double min=0;

            // find the minimum product of alpha, gamma, beta
            for (int u=0;u<2;u++)
                for (int state=0;state<numstates;state++)
                {
                    double temp=a[0][k][state]*gammaE[0][u][k][state]*b[0][k][tostate[u][state]];
                    if ((temp<min&&temp!=0)|| min==0)
                        min=temp;
                }

            // if all else fails, make min real small
            if (min == 0 || min > 1)
                min=1E-100;

            double topbottom[2];
            for (int u=0;u<2;u++)
            {
                topbottom[u]=0.0;
                for (int state=0;state<numstates;state++)
                    topbottom[u]+=(a[0][k][state]*gammaE[0][u][k][state]*b[0][k][tostate[u][state]]);
            }

            if (topbottom[0] == 0)
                topbottom[0]=min;
            else if (topbottom[1] == 0)
                topbottom[1]=min;

            L[0][k]=(log(topbottom[1]/topbottom[0]));
        }

        interleavedouble(L[0],size);
        // remember to deinterleave for next iteration
        interleavedouble(mesg,size);

        // start decoder 2
        // code almost same as decoder 1,
        // could combine code into one but too lazy
        for (int k=0;k<size;k++)
        {
            // calculate the gamma(s',s);
            for (int input=0;input<2;input++)
            {
                double uk=(input == 0) ? -1.0 : 1.0;

                for (int s=0;s<numstates;s++)
                {
                    double xk=(output[input][s] == 0) ? -1.0 : 1.0;

                    gammaE[1][input][k][s]=exp(0.5*Lc*parity2[k]*xk);
                    gamma[1][input][k][s]=exp(0.5*uk*(L[0][k]+Lc*mesg[k]))*gammaE[1][input][k][s];
                }
            }
        }

        // calculate the alpha terms
        for (int k=1;k<size;k++)
        {
            double temp=0;

            // calculate denominator
            for (int state=0;state<numstates;state++)
                temp += a[1][k-1][fromstate[0][state]]*gamma[1][0][k-1][fromstate[0][state]]+a[1][k-1][fromstate[1][state]]*gamma[1][1][k-1][fromstate[1][state]];

            for (int state=0;state<numstates;state++)
                a[1][k][state]=(a[1][k-1][fromstate[0][state]]*gamma[1][0][k-1][fromstate[0][state]]+a[1][k-1][fromstate[1][state]]*gamma[1][1][k-1][fromstate[1][state]])/temp;
        }

        // in the first iteration, set b[1][N-1]=a[1][N-1] for decoder 2
        // this decoder can't be terminated to state 0 because of the
        // interleaver. The performance loss is supposedly neglible.
        if (c==0)
        {
            double temp=0;

            // calculate denominator
            for (int state=0;state<numstates;state++)
                temp+=a[1][size-1][fromstate[0][state]]*gamma[1][0][size-1][fromstate[0][state]]+a[1][size-1][fromstate[1][state]]*gamma[1][1][size-1][fromstate[1][state]];

            for (int state=0;state<numstates;state++)
                b[1][size-1][state]=(a[1][size-1][fromstate[0][state]]*gamma[1][0][size-1][fromstate[0][state]]+a[1][size-1][fromstate[1][state]]*gamma[1][1][size-1][fromstate[1][state]])/temp;
        }

        for (int k=size-1;k>=1;k--)
        {
            double temp=0;

            // calculate denominator
            for (int state=0;state<numstates;state++)
                temp += a[1][k][fromstate[0][state]]*gamma[1][0][k][fromstate[0][state]]+a[1][k][fromstate[1][state]]*gamma[1][1][k][fromstate[1][state]];

            for (int state=0;state<numstates;state++)
                b[1][k-1][state]=(b[1][k][tostate[0][state]]*gamma[1][0][k][state]+b[1][k][tostate[1][state]]*gamma[1][1][k][state])/temp;
        }

        for (int k=0;k<size;k++)
        {
            double min=0;

            // find the minimum product of alpha, gamma, beta
            for (int u=0;u<2;u++)
                for (int state=0;state<numstates;state++)
                {
                    double temp=a[1][k][state]*gammaE[1][u][k][state]*b[1][k][tostate[u][state]];

                    if ((temp < min && temp != 0)|| min == 0)
                        min=temp;
                }
            // if all else fails, make min real small
            if (min == 0 || min > 1)
                min=1E-100;

            double topbottom[2];

            for (int u=0;u<2;u++)
            {
                topbottom[u]=0.0;

                for (int state=0;state<numstates;state++)
                    topbottom[u] += (a[1][k][state]*gammaE[1][u][k][state]*b[1][k][tostate[u][state]]);
            }

            if (topbottom[0] == 0)
                topbottom[0]=min;
            else if (topbottom[1] == 0)
                topbottom[1]=min;

            L[1][k]=(log(topbottom[1]/topbottom[0]));
        }

        deinterleavedouble(mesg,size);
        deinterleavedouble(L[1],size);
        // get L[0] back to normal after decoder 2
        deinterleavedouble(L[0],size);

        bool temp=true;
        for (int k=0;k<size;k++)
            if (boolmesg[k]!=((Lc*mesg[k]+L[0][k]+L[1][k]>0.0)?true:false))
                temp=false;

        // we can quit prematurely since it has been decoded
        if (temp==true)
        {
            returnvalue=c+1;
            c=ITERATIONS;
        }
    }
    // end decoder 2

    // make decisions
    for (int k=0;k<size;k++)
        if ((Lc*mesg[k] + L[0][k] + L[1][k]) > 0)
            mesg[k]=1.0;
        else
            mesg[k]=-1.0;

    return returnvalue;
}




void deinterleavedouble(double *mesg, unsigned size)
{
    double *temp;

    temp=new double[size];

    for (int x=0;x<size;x++)
        temp[x]=mesg[x];

    for (int x=0;x<size;x++)
        mesg[deinterleavearray[x]]=temp[x];

    delete [] temp;
}




void interleavedouble(double *mesg, unsigned size)
{
    double *temp;

    temp=new double[size];

    for (int x=0;x<size;x++)
        temp[x]=mesg[x];

    for (int x=0;x<size;x++)
        mesg[interleavearray[x]]=temp[x];

    delete [] temp;
}




void interleave(bool *mesg,unsigned size)
{
    bool *temp;

    temp=new bool[size];

    for (int x=0;x<size;x++)
        temp[x]=mesg[x];

    for (int x=0;x<size;x++)
        mesg[interleavearray[x]]=temp[x];

    delete [] temp;
}




void deinterleave(bool *mesg,unsigned size)
{
    bool *temp;

    temp=new bool[size];

    for (int x=0;x<size;x++)
        temp[x]=mesg[x];

    for (int x=0;x<size;x++)
        mesg[deinterleavearray[x]]=temp[x];

    delete [] temp;
}




void createinterleave(unsigned size)
{
    bool *yesno;

    yesno=new bool[size];

    for (int x=0;x<N;x++)
        yesno[x]=false;

    // create an interleave array
    for (int x=0;x<N;x++)
    {
        unsigned  val;

        do
        {
            val=before.longrandom(N);
        }
        while (yesno[val] == true);

        yesno[val]=true;
        interleavearray[x]=val;
        deinterleavearray[val]=x;
    }

    delete [] yesno;
}




void encode(bool *mesg,bool *parity,unsigned size, bool force)
{

    unsigned state=0;

    for (int x=0;x<size;x++)
    {
        // force the encoder to zero state at the end
        if (x>=size-memory && force)
        {
            if (tostate[0][state]&1)
                mesg[x]=true;
            else
                mesg[x]=false;
        }

        // can't assume the bool type has an intrinsic value of 0 or 1
        // may differ from platform to platform
        int uk=mesg[x] ?1:0;

        // calculate output due to new mesg bit
        parity[x]=output[uk][state];
        // calculate the new state
        state=tostate[uk][state];
    }
}


bool add(bool a, bool b)
{
    return a==b ? false:true;
}


double gaussian(double sigma)
{
    double rndm, u1, u2, s, noise;
    do
    {
        rndm=(double)rand()/(double)RAND_MAX;
        u1=rndm*2-1.0;
        rndm=(double)rand()/(double)RAND_MAX;
        u2=rndm*2-1.0;
        s=u1 * u1 + u2 * u2;
    }
    while ( s >= 1 );
    noise=u1 * sqrt( (-2.0*log(s))/s );
    return(noise*sigma);
}