turbo.cpp

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

#include <cstdio>
#include <cmath>
#include <cstring>
#include <cstdlib>
#include <climits>
#include "random.h"

int          states;
long         seed;
unsigned int N;
unsigned int ITERATIONS;
long         maxsim;
double       No;
double       sigma;
int          temp;

double gaussian(double variance);
void encode(bool *mesg,bool *parity1,unsigned size, bool force);
void interleave(bool *mesg, unsigned size);
bool add(bool a,bool b);
void interleave(bool *mesg, unsigned size);
void deinterleave(bool *mesg, unsigned size);
void deinterleavedouble(double *mesg, unsigned size);
void interleavedouble(double *mesg, unsigned size);
void createinterleave(unsigned size);
void createencodetable(int states);
void inttobool(unsigned state, bool *array, unsigned size);
void booltoint(bool *array, unsigned size, unsigned *state);
unsigned decode (double *channelmesg, double *parity1, double *parity2, unsigned size, bool *mesg);

unsigned *interleavearray;
unsigned *deinterleavearray;

unsigned numstates;
unsigned memory;

// [2]=input, [16]=current state, tostate[2][16]=next state
unsigned *tostate[2];

// [2]=last input, [16]=current state, fromstate[2][16]=previous state
unsigned *fromstate[2];

// [2]=input, [16]=current state, output[2][16]=output of encoder
bool *output[2];

Random before;

int main(int argc, char *argv[])
{
    bool *mesg;
    bool *parity1;
    bool *parity2;
    double *channelmesg;
    double *channelparity1;
    double *channelparity2;
    register int i;
    char name1[40], name2[40];
    FILE *fp;

    // Command line processing
    if (argc != 10)
    {
        printf("Usage: %s Num_states Max_iter INT_file", argv[0]);
        printf(" BER_file EbNo_init EbNo_final EbNo_inc  max_sims  seed\n");
        exit(0);
    }

    double INIT, FINAL, INC;

    sscanf(argv[1],"%d",  &states);
    sscanf(argv[2],"%d",  &ITERATIONS);
    sscanf(argv[3],"%s",  name1);
    sscanf(argv[4],"%s",  name2);
    sscanf(argv[5],"%lf", &INIT);
    sscanf(argv[6],"%lf", &FINAL);
    sscanf(argv[7],"%lf", &INC);
    sscanf(argv[8],"%ld", &maxsim);
    sscanf(argv[9],"%ld", &seed);



    // READ INTERLEAVER
    fp=fopen(name1,"r");
    N=0;
    while (!(fscanf(fp,"%d", &temp)==EOF))
        N++;
    fclose(fp);
    interleavearray=new unsigned[N];
    deinterleavearray=new unsigned[N];
    fp=fopen(name1,"r");
    i=0;
    while (!(fscanf(fp,"%d", &interleavearray[i])==EOF))
        i++;
    fclose(fp);
    for (i=0; i<N; i++) deinterleavearray[interleavearray[i]]=i;

    double rate=(double)N / (3.0*(double)N + 4.0);

    printf("\nParallel concatenated (turbo) code of rate=%lf\n", rate);
    printf("%2d-state constituent recursive convolutional codes \n\n", states);

    Random(seed);

    mesg=new bool[N];
    parity1=new bool[N];
    parity2=new bool[N];
    channelmesg=new double[N];
    channelparity1=new double[N];
    channelparity2=new double[N];


    fp=fopen(name2,"w");
    createencodetable(states);

    for (double dB=INIT; dB<FINAL+1.0e-1; dB += INC)
    {
        unsigned totalN=0;
        unsigned totalerrors=0;
        unsigned numiter;
        unsigned totaliter=0;
        unsigned numloops=0;

        // dB=10*log(Eb/No) where Eb is 1
        // No=1.0/pow(10.0,dB/10.0);
        sigma=1.0/sqrt(2.0*0.3333333*pow(10.0,dB/10.0));

        bool keepgoing=true;

        do
        {
            for (int x=0;x<N;x++)
                mesg[x]=before.boolrandom();

            // create a random interleaver for each decode trial
            // createinterleave(N);

            encode(mesg,parity1,N,true);
            interleave(mesg,N);
            encode(mesg,parity2,N,false);
            deinterleave(mesg,N);

            for (int x=0;x<N;x++)
            {
                mesg[x]    ?    channelmesg[x]=1.0 :    channelmesg[x]=-1.0;
                parity1[x] ? channelparity1[x]=1.0 : channelparity1[x]=-1.0;
                parity2[x] ? channelparity2[x]=1.0 : channelparity2[x]=-1.0;
            }

            // add gaussian noise, mean=0, variance=No/2
            for (int x=0;x<N;x++)
            {
                //channelmesg[x] += gaussian(No/2);
                //channelparity1[x] += gaussian(No/2);
                //channelparity2[x] += gaussian(No/2);
                channelmesg[x] += gaussian(sigma);
                channelparity1[x] += gaussian(sigma);
                channelparity2[x] += gaussian(sigma);
            }

            numiter=decode(channelmesg, channelparity1, channelparity2,N,mesg);

            unsigned numerrors=0;

            for (int x=0;x<N;x++)
            {
                bool temp=channelmesg[x] == 1 ? true : false;

                if (mesg[x] != temp)
                    numerrors++;
            }

            totalerrors += numerrors;
            totalN += N;
            totaliter += numiter;
            numloops++;

            char name[30];

            // stopping criterion
            if (totalerrors> 5000)
                keepgoing=false;
            else if (totalN>maxsim)
                keepgoing=false;
        }
        while (keepgoing == true);

        printf ("%5.2lf  %10.7e  %12ld %12ld     %lf     %d\n",dB,(totalerrors/(double)totalN),totalerrors,totalN,(double)totaliter/(double)numloops, N);
        fprintf(fp, "%8.4lf %10.7e\n", dB, (totalerrors/(double)totalN));
        fflush(fp);
    }

    delete [] interleavearray;
    interleavearray=NULL;
    delete [] deinterleavearray;
    deinterleavearray=NULL;
    delete [] mesg;
    mesg=NULL;
    delete [] parity1;
    parity1=NULL;
    delete [] parity2;
    parity2=NULL;
    delete [] channelmesg;
    channelmesg=NULL;
    delete [] channelparity1;
    channelparity1=NULL;
    delete [] channelparity2;
    channelparity2=NULL;
}

void createencodetable(int states)
{
    bool *boolstate;
    bool *newstate;

    if (states == 4)           // Generators {7,5}
    {
        numstates=4;
        memory=2;

        // create arrays used by encode and decode
        output[0]=new bool[numstates];
        output[1]=new bool[numstates];
        fromstate[0]=new unsigned[numstates];
        fromstate[1]=new unsigned[numstates];
        tostate[0]=new unsigned[numstates];
        tostate[1]=new unsigned[numstates];

        boolstate=new bool[memory];
        newstate=new bool[memory];

        for (unsigned input=0;input<2;input++)
            for (unsigned intstate=0;intstate<numstates;intstate++)
            {
                bool boolinput=(input==0)?false:true;

                inttobool(intstate,boolstate,memory);

                // calculate output due to the input
                output[input][intstate]=add(boolinput,boolstate[0]);
                output[input][intstate]=add(output[input][intstate],boolstate[1]);

                // calculate new states
                newstate[1]=boolstate[0];
                newstate[0]=add(add(boolinput,boolstate[0]),boolstate[1]);
                // from s' to s
                booltoint (newstate,memory,&tostate[input][intstate]);
                // from s to s'
                fromstate[input][tostate[input][intstate]]=intstate;
            }

    }

    else if (states == 8)           // Generators {17,15}
    {
        numstates=8;
        memory=3;

        // create arrays used by encode and decode
        output[0]=new bool[numstates];
        output[1]=new bool[numstates];
        fromstate[0]=new unsigned[numstates];
        fromstate[1]=new unsigned[numstates];
        tostate[0]=new unsigned[numstates];
        tostate[1]=new unsigned[numstates];

        boolstate=new bool[memory];
        newstate=new bool[memory];

        for (unsigned input=0;input<2;input++)
            for (unsigned intstate=0;intstate<numstates;intstate++)
            {
                bool boolinput=(input == 0) ?false:true;

                inttobool(intstate,boolstate,memory);

                // calculate output due to the input
                //output[input][intstate]=add(boolinput,boolstate[0]);
                output[input][intstate]=add(boolinput,boolstate[1]);
                //output[input][intstate]=add(output[input][intstate],boolstate[2]);
                output[input][intstate]=add(output[input][intstate],boolstate[2]);

                // calculate new states
                newstate[2]=boolstate[1];
                newstate[1]=boolstate[0];
                newstate[0]=add(add(add(boolinput,boolstate[0]),boolstate[1]),boolstate[2]);
                // from s' to s
                booltoint (newstate,memory,&tostate[input][intstate]);
                // from s to s'
                fromstate[input][tostate[input][intstate]]=intstate;
            }

    }

    else if (states == 16)           // Generators {31,27}
    {
        numstates=16;
        memory=4;

        // create arrays used by encode and decode
        output[0]=new bool[numstates];
        output[1]=new bool[numstates];
        fromstate[0]=new unsigned[numstates];
        fromstate[1]=new unsigned[numstates];
        tostate[0]=new unsigned[numstates];
        tostate[1]=new unsigned[numstates];

        boolstate=new bool[memory];
        newstate=new bool[memory];

        for (unsigned input=0;input<2;input++)
            for (unsigned intstate=0;intstate<numstates;intstate++)
            {
                bool boolinput=(input == 0) ? false : true;

                inttobool(intstate,boolstate,memory);

                // calculate output due to the input
                output[input][intstate]=add(boolinput,boolstate[0]);
                output[input][intstate]=add(output[input][intstate],boolstate[3]);
                output[input][intstate]=add(output[input][intstate],boolstate[1]);
                output[input][intstate]=add(output[input][intstate],boolstate[2]);
                output[input][intstate]=add(output[input][intstate],boolstate[3]);

                // calculate new states
                newstate[3]=boolstate[2];
                newstate[2]=boolstate[1];
                newstate[1]=boolstate[0];
                newstate[0]=add(add(boolinput,boolstate[0]),boolstate[3]);
                // from s' to s
                booltoint (newstate,memory,&tostate[input][intstate]);
                // from s to s'
                fromstate[input][tostate[input][intstate]]=intstate;
            }

    }

    delete [] boolstate;
    delete [] newstate;
}





void inttobool(unsigned state, bool *array, unsigned size)
{
    for (unsigned x=0;x<size;x++)
    {
        unsigned next=state >> 1;

        if ((next << 1) == state)
            array[x]=false;
        else
            array[x]=true;
        state=next;
    }
}


void booltoint(bool *array, unsigned size, unsigned *state)
{
    *state=0;
    for (int x=0;x<size;x++)
        if (array[x] == true)
            (*state) |= (1 << x);
}


unsigned decode (double *mesg, double *parity1, double *parity2, unsigned size, bool *boolmesg)
{
    static double **a[2];
    static double **b[2];
    static double *L[2];
    static double **gamma[2][2];
    static double **gammaE[2][2];
    static bool initialized=false;
    static unsigned lastsize=0;
    unsigned returnvalue=ITERATIONS;

    // minimize new's and delete's to only when needed
    if (size != lastsize && initialized == true)
    {
        // delete all the arrays and rebuild
        for (int y=0;y<2;y++)
            for (int x=0;x<2;x++)
            {
                for (int z=0;z<lastsize;z++)
                {
                    delete [] gamma[y][x][z];
                    gamma[y][x][z]=NULL;
                    delete [] gammaE[y][x][z];
                    gammaE[y][x][z]=NULL;
                }

                delete [] gamma[y][x];
                gamma[y][x]=NULL;
                delete [] gammaE[y][x];
                gammaE[y][x]=NULL;
            }

        // create L[encoder #]
        for (int y=0;y<2;y++)
        {
            delete [] L[y];
            L[y]=NULL;
        }
        // create alpha[encoder #][k][state]
        for (int x=0;x<2;x++)
        {
            for (int y=0;y<lastsize;y++)
            {
                delete [] a[x][y];
                a[x][y]=NULL;
                delete [] b[x][y];
                b[x][y]=NULL;
            }

            delete [] a[x];
            a[x]=NULL;
            delete [] b[x];