soft_viterbi.c

瑞利衰落信道的Matlab仿真模型程序

/*=================================================================
 *
 * soft_viterbi.C	.MEX file corresponding to soft_viterbi.M
 *	        
 * The calling syntax is:
 *
 *		 [out_put] = soft_viterbi(channel_output)
 *
 *  You may also want to look at the corresponding M-code, soft_viterbi.m.
 *
 * This is a MEX-file for MATLAB.  
 * Copyright 2005 The MathWorks, Inc.
 *
 *=================================================================*/

#include "mex.h"
#include "matrix.h"
#include <memory.h>
#include <stdlib.h>
#include <malloc.h>
#include <stdio.h>

#define	ChannelOutput_IN	prhs[0]

#define	Output_OUT	plhs[0]


void soft_viterbi(double channel_output[], double out_put[])

{

	int result[4992] = {0};
	const num_of_states = 256;
	int dec_1,dec_2;
	const depth_of_trellis = 4992;						//depth_of_trellis:需要译码的次数
	int *rout = (int *)malloc(5111808);					//save path定义存储路径变量的数组,长度为depth_of_trellis
	int branch_output_1[2];							//对应不同的寄存器值的信道编码后的输出码组；
	int branch_output_2[2];							//对应不同的寄存器值的信道编码后的输出码组；
	int c=0;
	int G[18] = {1,0,0,0,1,1,1,0,1,1,1,0,1,0,1,1,1,1};			//生成矩阵G
	const L = 9;
	int i,j,q;  
	int x,y;
	double temp1 = 0;
	double temp2 = 0;
	int d=0;										
	int sum0=0;
	int sum1=0;
	int sum2=0;
	int sum3=0;
	int register_1[2304]={0};					//num_of_states*L = 2304
	int register_2[2304]={0};
	double tolerance_1[256]={0};				
	double tolerance_2[256]={0};
	int current_state = 0;



	//进行初始化工作；
	//对寄存器register_1、register_2赋值；
	for	(i = 0;i < num_of_states;i++)
	{
		dec_1 = num_of_states + i;
		dec_2 = i;
		//第2至第L列初始化为状态数的二进制形式；
		for (j = L-1 ;j >= 0;j--)
		{
			register_1[i*9+j] = dec_1%2;
			dec_1 = (int)(dec_1/2);
			register_2[i*9+j] = dec_2%2;
			dec_2 = (int)(dec_2/2);
		}
		tolerance_1[i] = -64;
	}
	tolerance_1[0] = 0;


	/*fp = fopen( "e:\\data","rt");
    if (fp != NULL)
	{
		aa = fread(channel_output,sizeof(int),9984,fp);
	}*/

	/*srand( (unsigned)time( NULL ) );	
	for( c = 0; c < 9984;c++ )
	{
		 channel_output[c] = rand() ;
		 if	( channel_output[c] > 16383 )
		 {
			 channel_output[c] = 2;
		 }
		 else
		 {
			 channel_output[c] = -2;
		 }

	}
	fp = fopen("myfirst.txt","w");
	if (fp != NULL)
	{
		for (x=0;x<9984;x++)
		{
			fprintf(fp,"%f ",channel_output[x]);
		}
	}
	
	fclose(fp);*/

	//fp=NULL;


	//对channel_output进行循环译码；
	for (j = 0;j<depth_of_trellis;j++)
	{
		//对register_2所有状态+输入求其对数似然比之和;
		for (i = 0;i<num_of_states;i++)                                          
		{                                                                        
			                                                    
			temp2 = 0;                                                       
			for(d=0;d<L;d++)                                                 
			{                                                                
				sum0=register_2[i*9+d]*G[d]+sum0;                      
				sum1=register_2[i*9+d]*G[9+d]+sum1;                     
			}            
			//branch_output:对应不同的寄存器值的信道编码后的输出码组；                                                    
			branch_output_2[0]=sum0%2;	                                 
			branch_output_2[1]=sum1%2;                                      
                                                                                         
			sum0=0;                                                          
			sum1=0;                                                          
                                                                                         
		//对似然比进行累加；                                                                
        	for(q=0;q<2;q++)                                                         
			{                   
				//如果为0，给似然比赋正值；                                             
				if(branch_output_2[q] == 0)                              
				{                                                        
					temp2 = temp2 + channel_output[j*2+q];           
				}    
				//如果为1，给似然比赋负值；                                                    
				else                                                     
				{                                                        
					temp2 = temp2 - channel_output[j*2+q];           
				}                                              
				                                                         
			}                                                                
			      
			//register_2中的度量值来自于register_1中的相应度量值；
			tolerance_2[i] = tolerance_1[(int)(i/2)] + temp2;         
		}    	
		
		
		//对register_1所有状态+输入求其对数似然比之和；
		for (i = 0;i<num_of_states;i++)                                          
		{                                                                        
			                                                                 
			//temp1 = 0;                                                       
			for(d=0;d<L;d++)                                                 
			{                                                             
				sum2=register_1[i*9+d]*G[d]+sum2;                      
				sum3=register_1[i*9+d]*G[9+d]+sum3;                      
			}       
			//branch_output:对应不同的寄存器值的信道编码后的输出码组；                                                         
			branch_output_1[0]=sum2%2;                                       
			branch_output_1[1]=sum3%2;                                       
                                                                                        
			sum2=0;                                                          
			sum3=0;                                                          
                                                                                         
			      
		//对似然比进行累加；
        	for(q=0;q<2;q++)                                                         
			{            
				//如果为0，给似然比赋正值；                                             
				if(branch_output_1[q] == 0)                             
				{                                                        
					temp1 = temp1 + channel_output[j*2+q];           
				}    
				//如果为1，给似然比赋负值；                                                    
				else                                                     
				{                                                        
					temp1 = temp1 - channel_output[j*2+q];           
				}                                                        
				                                                         
			}                                                                
		   	//register_1中的度量值来自于register_1中的相应度量值；
			tolerance_1[i] = tolerance_1[(int)(i/2)+num_of_states/2] + temp1;
            temp1 = 0;                                                                              
		}                                                                                                                                            

		//对度量值进行判断，如果register_1中的度量值较小，则将其更改为register_2中的度量值；
		for(i=0;i<num_of_states;i++)
		{
		
			if(tolerance_1[i] < tolerance_2[i])  
			{
				tolerance_1[i] = tolerance_2[i];
				//同时将路径赋值rout中的内容也更新；
				rout[i*4992+j] = register_2[i*9];
			}
			else									
			{
				//记录相应的register_1中的路径信息；
				rout[i*4992+j] = register_1[i*9];	
			}
		}			
		
	}


	//trace back(回溯部分)


	//对深度为depth_of_trellis的rout中的内容进行回溯；
	for(j=(depth_of_trellis-1);j>=0;j--)
	{
		//result记录回溯过程中找出的路径信息；
		result[j] = current_state%2;   
			
		//如rout(current_state,j)为0，说明下一状态值在所有状态的前半部分；
		if(rout[current_state*4992+j] == 0)		//下一状态值为当前状态值的一半；
			current_state = (int)(current_state/2);
		//下一状态值为当前状态值的一半加上128；
		else
			current_state = (int)(current_state/2)+num_of_states/2;	

	}
	//输出结果；
	for(i=0;i<4984;i++)
	{
		out_put[i]=(double)(result[i]);					//?
	}
	
	
	/*fp_2 = fopen("mysecond.txt","w");
	if (fp_2 != NULL)
	{
		for (y=0;y<4984;y++)
		{
			fprintf(fp_2,"%f\n",out_put[y]);
		}
	}
	fclose(fp_2);*/
	free(rout);
	return ;
}
//mex文件接口部分程序
//mexFunction
void mexFunction( int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[] )
{ 
  double *out_put; 
  double *channel_output; 
    
    /*Check for proper number of arguments */
    
    if (nrhs != 1) 
	{
		mexErrMsgTxt("One input arguments required."); 
    }
	else if (nlhs > 1) 
	{
		mexErrMsgTxt("Too many output arguments."); 
	} 
    
    
    /* Create a matrix for the return argument */
   Output_OUT = mxCreateDoubleMatrix(1, 4984, mxREAL); 
    
    /* Assign pointers to the various parameters */
   out_put = mxGetPr(Output_OUT);
    
   channel_output = mxGetPr(ChannelOutput_IN); 
        
    /* Do the actual computations in a subroutine */
   soft_viterbi(channel_output,out_put); 
    
}