decode.c

这是一个vblast程序

#include <stdio.h>
#include <stdlib.h>
#include "matlab.h"
#include "mex.h"
#include "math.h"
#include "matrix.h"
#define pi 3.1415926
#define block_size 800
#define no_of_states 8

mxArray *BCJR(mxArray **LLR_I,mxArray *In_LLR,mxArray *code_g)
{   
   double alpha[block_size+1][no_of_states],branch[no_of_states],alpha_sum[block_size+1]={0},belta[block_size+1][no_of_states];
   double den[no_of_states],num[no_of_states];
   double den_c1[no_of_states],num_c1[no_of_states];
   double den_c2[no_of_states],num_c2[no_of_states];
   double *Data_p,*nxt_s_p,*nxt_o_p;
   double branch0,branch1,den_sum,num_sum,sum1,llr[block_size];
   double branch0_c1,branch1_c1,den_sum_c1,num_sum_c1,sum1_c1,llr_c1[block_size],temp0_c1,temp1_c1;
   double branch0_c2,branch1_c2,den_sum_c2,num_sum_c2,sum1_c2,llr_c2[block_size],temp0_c2,temp1_c2;
   double temp_c1,temp_c2;
   int i1,cnt_s,j1,j2;
   mxArray *lhs[2],*LLR1_C1=NULL,*LLR1_C2=NULL,*LLR_C=NULL;
   lhs[0]=NULL;
   lhs[1]=NULL;

   mlfEnterNewContext(1,2,LLR_I,In_LLR,code_g);
   
   Data_p=mxGetPr(In_LLR);

   
   mexCallMATLAB(2,lhs,1,&code_g,"cnc_trellis"); 
   
   nxt_s_p=mxGetPr(lhs[0]);
   nxt_o_p=mxGetPr(lhs[1]);

   alpha[0][0]=1;
   for(j1=1;j1<no_of_states;j1++)
   {
           alpha[0][j1]=0; 
   }
   
   /* forward recursion */
   
   for(i1=0;i1<block_size;i1++)
     {
        for(cnt_s=0;cnt_s<no_of_states;cnt_s++) 
          {
              for(j1=0;j1<no_of_states;j1++)
              {
                 branch[j1]=0;
               }
               
              /* input bit */ 
              
              for(j2=0;j2<no_of_states;j2++) 
               {    
                 temp_c1=0;
                 temp_c2=0;
                 if(cnt_s==(int)(*(nxt_s_p+j2)))         /* input 0 */
                   {
                 
                     /* output c1 */
                 
                     if(*(nxt_o_p+j2)==0)
                      {
                         temp_c1=1/(1+exp((*(Data_p+2*i1))));
                       }
                     else
                      { 
                         temp_c1=exp((*(Data_p+2*i1)))/(1+exp((*(Data_p+2*i1))));
                       }
                      
                     /* output c2 */  
                       
                     if(*(nxt_o_p+no_of_states+j2)==0)
                      {
                         temp_c2=1/(1+exp((*(Data_p+2*i1+1))));
                       }
                     else
                      { 
                         temp_c2=exp((*(Data_p+2*i1+1)))/(1+exp((*(Data_p+2*i1+1))));
                       }
                     branch[j2]=branch[j2]+temp_c1*temp_c2;
                   }
                  
                  if(cnt_s==(int)(*(nxt_s_p+no_of_states+j2)))         /* input 1 */
                    {
                 
                     /* output c1 */
                 
                     if(*(nxt_o_p+2*no_of_states+j2)==0)
                      {
                         temp_c1=1/(1+exp((*(Data_p+2*i1))));
                       }
                     else
                      { 
                         temp_c1=exp((*(Data_p+2*i1)))/(1+exp((*(Data_p+2*i1))));
                       }
                      
                     /* output c2 */  
                       
                     if(*(nxt_o_p+3*no_of_states+j2)==0)
                      {
                         temp_c2=1/(1+exp((*(Data_p+2*i1+1))));
                       }
                     else
                      { 
                         temp_c2=exp((*(Data_p+2*i1+1)))/(1+exp((*(Data_p+2*i1+1))));
                       }
                     branch[j2]=branch[j2]+temp_c1*temp_c2;
                   }
               }    
              
              sum1=0;
              for(j1=0;j1<no_of_states;j1++)
              {
                 sum1=sum1+branch[j1]*alpha[i1][j1];
               }
              if(sum1!=0)
                 alpha[i1+1][cnt_s]=sum1;  
              else alpha[i1+1][cnt_s]=0;
           }
        for(j1=0;j1<no_of_states;j1++)
        {
           alpha_sum[i1+1]=alpha_sum[i1+1]+alpha[i1+1][j1];
         }
        for(j1=0;j1<no_of_states;j1++)
        {
           alpha[i1+1][j1]=alpha[i1+1][j1]/alpha_sum[i1+1];
         }
      }     


      
   /* backward recursion */
   

    belta[block_size-1][0]=1;

    for(j1=1;j1<no_of_states;j1++)
        {
           belta[block_size-1][j1]=0;
        }


   
   for(i1=block_size-2;i1>=0;i1--)
     {
       for(cnt_s=0;cnt_s<no_of_states;cnt_s++)
         {
            for(j1=0;j1<no_of_states;j1++)
            {
                 branch[j1]=0;
             }

             /* input bit 0 */
             
             /* output c1 */
             
            temp_c1=0;
            temp_c2=0;
            if(*(nxt_o_p+cnt_s)==0)
              {
                temp_c1=1/(1+exp((*(Data_p+2*(i1+1)))));
               }
            else
              {
                 temp_c1=exp((*(Data_p+2*(i1+1))))/(1+exp((*(Data_p+2*(i1+1)))));
               }
               
            /* output c2 */   
            
            if(*(nxt_o_p+no_of_states+cnt_s)==0)
              {
                temp_c2=1/(1+exp((*(Data_p+2*(i1+1)+1))));
               }
            else
              {
                 temp_c2=exp((*(Data_p+2*(i1+1)+1)))/(1+exp((*(Data_p+2*(i1+1)+1))));
               } 
               
            branch[(int)(*(nxt_s_p+cnt_s))]=temp_c1*temp_c2;  
            
             /* input bit 1 */
             
             /* output c1 */
             
             temp_c1=0;
             temp_c2=0;
            if(*(nxt_o_p+2*no_of_states+cnt_s)==0)
              {
                temp_c1=1/(1+exp((*(Data_p+2*(i1+1)))));
               }
            else
              {
                 temp_c1=exp((*(Data_p+2*(i1+1))))/(1+exp((*(Data_p+2*(i1+1)))));
               }
               
            /* output c2 */  
            
            if(*(nxt_o_p+3*no_of_states+cnt_s)==0)
              {
                temp_c2=1/(1+exp((*(Data_p+2*(i1+1)+1))));
               }
            else
              {
                 temp_c2=exp((*(Data_p+2*(i1+1)+1)))/(1+exp((*(Data_p+2*(i1+1)+1))));
               } 
               
            branch[(int)(*(nxt_s_p+no_of_states+cnt_s))]=branch[(int)(*(nxt_s_p+no_of_states+cnt_s))]+temp_c1*temp_c2;

            sum1=0;
            for(j1=0;j1<no_of_states;j1++)
            {
               sum1=sum1+branch[j1]*belta[i1+1][j1];   
             }
            if(sum1!=0)
               belta[i1][cnt_s]=sum1;  
            else belta[i1][cnt_s]=0;
          }      
        for(j1=0;j1<no_of_states;j1++)
        {
           belta[i1][j1]=belta[i1][j1]/alpha_sum[i1+1];
        }
     }
     
     
     
  for(i1=0;i1<block_size;i1++)
   {
     for(cnt_s=0;cnt_s<no_of_states;cnt_s++)
     {
        den[cnt_s]=0;
        num[cnt_s]=0;
        den_c1[cnt_s]=0;
        num_c1[cnt_s]=0;
        den_c2[cnt_s]=0;
        num_c2[cnt_s]=0;
     }
     for(cnt_s=0;cnt_s<no_of_states;cnt_s++)
       {
       
          /* input bit */
          
           for(j2=0;j2<no_of_states;j2++) 
             {    
                branch0=0;
                branch1=0;
                 
                if(cnt_s==(int)(*(nxt_s_p+j2)))         /* input 0 */
                   {
                 
                     /* output c1 */
                 
                     if(*(nxt_o_p+j2)==0)
                      {
                         temp_c1=1/(1+exp((*(Data_p+2*i1))));
                       }
                     else
                      { 
                         temp_c1=exp((*(Data_p+2*i1)))/(1+exp((*(Data_p+2*i1))));
                       }
                      
                     /* output c2 */  
                       
                     if(*(nxt_o_p+no_of_states+j2)==0)
                      {
                         temp_c2=1/(1+exp((*(Data_p+2*i1+1))));
                       }
                     else
                      { 
                         temp_c2=exp((*(Data_p+2*i1+1)))/(1+exp((*(Data_p+2*i1+1))));