logmap.c

与之前一个对应的代码

#include<stdio.h>
#include<math.h>
#include"ctc.h"
void logmap(float rec[],double L_a_1row[],int iter,int puncture,int L1,double L_all[],int *start_state)
{
	double Alpha[4][L_total+1][nstates];
	double Beta[4][L_total+1][nstates];
	double tempmax[L_total+1];
	double gamma[4][nstates];
	double p_endstate[nstates];
	double p_startstate[nstates];
	double p2[nstates],p1;
	double temp[4][nstates];
	double delta[4];
	double L_a[2][L_total];
	int i,j,k,n;
	//initialization
	if(iter<L1){
		for(i=0;i<4;i++){
			for(j=0;j<nstates;j++){
				Alpha[i][0][j]=-log(nstates);
				Beta[i][L_total][j]=-log(nstates);
			}
		}
	}
	else{
		for(i=0;i<4;i++){
			for(j=0;j<nstates;j++){
				Alpha[i][0][j]=-Infty;
				Beta[i][L_total][j]=-Infty;
			}
			Alpha[i][0][*start_state]=0;
			Beta[i][L_total][*start_state]=0;
		}
	}
	for(i=0;i<L_total;i++){
		L_a[0][i]=L_a_1row[2*i];
		L_a[1][i]=L_a_1row[2*i+1];
	}
	// Trace forward, compute Alpha
	for(k=1;k<L_total+1;k++){
		for(i=0;i<4;i++){
			for(j=0;j<nstates;j++){
				gamma[i][j]=-Infty;
			}
		}
		for(i=0;i<nstates;i++){
			gamma[0][i] = -rec[3*k-3]-rec[3*k-2]+rec[3*k-1]*state_map[i][10]-log(1+exp(L_a[0][k-1]))-log(1+exp(L_a[1][k-1]));
			gamma[1][i] = -rec[3*k-3]+rec[3*k-2]+rec[3*k-1]*state_map[i][13]-log(1+exp(L_a[0][k-1]))-log(1+exp(L_a[1][k-1]))+L_a[1][k-1];
			gamma[2][i] = rec[3*k-3]-rec[3*k-2]+rec[3*k-1]*state_map[i][16]-log(1+exp(L_a[0][k-1]))-log(1+exp(L_a[1][k-1]))+L_a[0][k-1];
			gamma[3][i] = rec[3*k-3]+rec[3*k-2]+rec[3*k-1]*state_map[i][19]-log(1+exp(L_a[0][k-1]))-log(1+exp(L_a[1][k-1]))+L_a[1][k-1]+L_a[0][k-1];
		}
		for(i=0;i<4;i++){
			for(j=0;j<nstates;j++){
				Alpha[i][k][j]=exp(Alpha[0][k-1][state_map[j][i]]+gamma[i][j])+exp(Alpha[1][k-1][state_map[j][i]]+gamma[i][j])
                +exp(Alpha[2][k-1][state_map[j][i]]+gamma[i][j])+exp(Alpha[3][k-1][state_map[j][i]]+gamma[i][j]);
				if(Alpha[i][k][j]<1e-300){
					Alpha[i][k][j] = -Infty;
				}
				else{
					Alpha[i][k][j] = log(Alpha[i][k][j]);
				}
			}
		}
		tempmax[k]=Alpha[0][k][0];
		for(i=0;i<4;i++){
			for(j=0;j<nstates;j++){
				if(tempmax[k]<Alpha[i][k][j]){
					tempmax[k]=Alpha[i][k][j];
				}
			}
		}
		for(i=0;i<4;i++){
			for(j=0;j<nstates;j++){
				Alpha[i][k][j]-=tempmax[k];
			}
		}
	}
	// Trace backward, compute Beta
	for(k=L_total-1;k>=0;k--){
		for(i=0;i<4;i++){
			for(j=0;j<nstates;j++){
				gamma[i][j]=-Infty;
			}
		}
		for(i=0;i<nstates;i++){
			gamma[0][i] = -rec[3*k]-rec[3*k+1]+rec[3*k+2]*state_map[i][22]-log(1+exp(L_a[0][k]))-log(1+exp(L_a[1][k]));
			gamma[1][i] = -rec[3*k]+rec[3*k+1]+rec[3*k+2]*state_map[i][25]-log(1+exp(L_a[0][k]))-log(1+exp(L_a[1][k]))+L_a[1][k];
			gamma[2][i] = rec[3*k]-rec[3*k+1]+rec[3*k+2]*state_map[i][28]-log(1+exp(L_a[0][k]))-log(1+exp(L_a[1][k]))+L_a[0][k];
			gamma[3][i] = rec[3*k]+rec[3*k+1]+rec[3*k+2]*state_map[i][31]-log(1+exp(L_a[0][k]))-log(1+exp(L_a[1][k]))+L_a[1][k]+L_a[0][k];
		}
		for(i=0;i<4;i++){
			for(j=0;j<nstates;j++){
				Beta[i][k][j]=exp(Beta[0][k+1][state_map[j][i+4]]+gamma[i][j])+exp(Beta[1][k+1][state_map[j][i+4]]+gamma[i][j])
                +exp(Beta[2][k+1][state_map[j][i+4]]+gamma[i][j])+exp(Beta[3][k+1][state_map[j][i+4]]+gamma[i][j]);
				if(Beta[i][k][j]<1e-300){
					Beta[i][k][j] = -Infty;
				}
				else{
					Beta[i][k][j] = log(Beta[i][k][j]);
				}
			}
		}
		for(i=0;i<4;i++){
			for(j=0;j<nstates;j++){
				Beta[i][k][j]-=tempmax[k+1];
			}
		}
	}
//	printf("%f\n",Alpha[0][L_total][0]);
	// Compute the soft output, log-likelihood ratio of symbols in the frame
	for(k=0;k<L_total;k++){
		for(i=0;i<nstates;i++){
			delta[0] = -rec[3*k]-rec[3*k+1]+rec[3*k+2]*state_map[i][10]-log(1+exp(L_a[0][k]))-log(1+exp(L_a[1][k]));
			delta[1] = -rec[3*k]+rec[3*k+1]+rec[3*k+2]*state_map[i][13]-log(1+exp(L_a[0][k]))-log(1+exp(L_a[1][k]))+L_a[1][k];
			delta[2] = rec[3*k]-rec[3*k+1]+rec[3*k+2]*state_map[i][16]-log(1+exp(L_a[0][k]))-log(1+exp(L_a[1][k]))+L_a[0][k];
			delta[3] = rec[3*k]+rec[3*k+1]+rec[3*k+2]*state_map[i][19]-log(1+exp(L_a[0][k]))-log(1+exp(L_a[1][k]))+L_a[1][k]+L_a[0][k];
			temp[0][i]=0;
			temp[1][i]=0;
			temp[2][i]=0;
			temp[3][i]=0;
			for(j=0;j<4;j++){
				for(n=0;n<4;n++){
					temp[0][i]+=exp(delta[0]+Alpha[n][k][state_map[i][0]]+Beta[j][k+1][i]);
					temp[1][i]+=exp(delta[1]+Alpha[n][k][state_map[i][1]]+Beta[j][k+1][i]);
					temp[2][i]+=exp(delta[2]+Alpha[n][k][state_map[i][2]]+Beta[j][k+1][i]);
					temp[3][i]+=exp(delta[3]+Alpha[n][k][state_map[i][3]]+Beta[j][k+1][i]);
				}
			}
		}
		delta[0]=temp[0][0];
		delta[1]=temp[1][0];
		delta[2]=temp[2][0];
		delta[3]=temp[3][0];
		for(i=1;i<nstates;i++){
			delta[0]+=temp[0][i];
			delta[1]+=temp[1][i];
			delta[2]+=temp[2][i];
			delta[3]+=temp[3][i];
		}
		L_all[2*k]=log(delta[2]+delta[3])-log(delta[0]+delta[1]);
		L_all[2*k+1]=log(delta[1]+delta[3])-log(delta[0]+delta[2]);
	}
	// estimate the start and end state
	if(iter==L1-1){
		p1=0;
		for(i=0;i<nstates;i++){
			p2[i]=0;
			for(j=0;j<4;j++){
				p2[i]+=exp(Alpha[j][L_total][i]);
			}
			p1+=p2[i];
		}
		for(i=0;i<nstates;i++){
			p_endstate[i]=log(p2[i])-log(p1);
		}
		p1=0;
		for(i=0;i<nstates;i++){
			p2[i]=0;
			for(j=0;j<4;j++){
				p2[i]+=exp(Beta[j][0][i]);
			}
			p1+=p2[i];
		}
		for(i=0;i<nstates;i++){
			p_startstate[i]=log(p2[i])-log(p1);
		}
		*start_state=0;
		for(i=0;i<nstates;i++){
			p2[i]=p_endstate[i]+p_startstate[i];
			if(p2[i]>p2[*start_state]){
				*start_state=i;
			}
		}
	}
}