polar SC译码(非递归) C语言+matlab混和编程

#include
#include
#include
#include "mex.h"
//声明函数
float f(float a,float b);
float g(float a,float b,int c);
int *decode(float *llr,int *CBR,int N,int n);
int mod(int a,int b);

//译码函数
int *decode(float *llr,int *CBR,int N,int n){
    //赋值llr CBR
    //mod执行次数
    int A[N]={};
    for(int i=0;i
        int j = i + 1;
        int layer = 0;
        while (j % 2 == 0) {
            j /= 2;
            layer += 1;
        }
        A[i] = layer;
    }
    //f执行次数
    int B[N] = {};
    B[0] = 0;
    for (int i = 1; i < N; i++) {
        int j = i;
        int layer = 0;
        while (j % 2 == 0) {
            j /= 2;
            layer += 1;
        }
        B[i] = layer;
    }
     //P数组------C数组
    double P[N][n + 1];
    int C[N][n + 1];
    //先把C[0][0]算出来，它比较特殊
    for (int i = 0; i < N; i++) {
        P[i][n] = llr[i];//赋值
    }
    int h = N;
    for (int j = n - 1; j >= 0; j--) {
        for (int i = 0; i < h / 2; i++) {
            P[i][j] = f(P[i][j + 1], P[i + h / 2][j + 1]);
        }
        h /= 2;
    }
if (!CBR[0]) { C[0][0] = 0; }
else { C[0][0] = (int)P[0][0] > 0 ? 0 : 1; }
//第二层开始
    for (int k = 1; k < N; k++) {
        if (k % 2 == 1) {//奇数执行g+mod
            if (!CBR[k]) { C[k][0] = 0; }
            else { float temp=(g(P[k - 1][1], P[k][1], C[k - 1][0]) > 0 ? 0 : 1);
                  C[k][0] =(int)temp;}
            //执行范围num1表示mod执行次数
            int num1 = A[k];
            //printf("n%d n",num1);
            //int flag = k;
            for (int j = 0; j 
                    // num2表示执行mod距离
                    int num2 = pow(2, j);
                    for (int l = 0; l < num2; l++) {
                        C[k - l - num2][j + 1] = mod(C[k - l - num2][j] ,C[k - l][j]);
                        C[k- l][j + 1] = C[k - l][j];
                    }
            }
        }
     if (k % 2 == 0) {//偶数执行g+f
            //num3表示f执行次数
            int num3 = B[k];
            //num4表示执行距离in
            //g运算
            int j = num3;
            int num4 = pow(2, num3);
            for (int r = 0; r < num4; r++) {
                P[k + r][j] = g(P[k - num4 + r][j + 1], P[k + r][j + 1], C[k - num4 + r][j]);
            }
            //下面开始执行f运算
			for (int q = j; q > 0; q--) {
                int num5 = pow(2, q - 1);
                for(int i=0;i
                    P[k+i][q - 1] = f(P[k+i][q], P[k + num5+i][q]);//害，一开始没有只写了k一层，P数组的相邻也要求一下！
                }
            }
            if (!CBR[k]) { C[k][0] = 0; }
            else
            {
                 float temp=P[k][0] > 0 ? 0 : 1;
                C[k][0]=(int)temp;
            }
        }
}
//返回译码值！
int *u_d;
u_d=(int *)malloc(sizeof(int)*N);
for(int i=0;i
    u_d[i]=C[i][0];
}
return u_d;
}


//f
float f(float a,float b) {
    float z = a * b > 0 ? 1 : -1;
    z = fmin(fabs(a), fabs(b)) * z;
    return z;
}
//g函数
float g(float a,float b,int c) {
    return (1 - 2 * c) * a + b;
}
//mod函数
int mod(int a,int b){
    if(a==b){
        return 0;
    }
    else{
        return 1;
    }
}
//接口函数
void mexFunction (int nlhs,mxArray *plhs[] ,int nrhs,const mxArray *prhs[]){
	//与matlab中function类似
	//右边输入需要四个参数、分别是11r、CBR(冻结比特)、N、n(log2(N))
    double *a=mxGetPr(prhs[0]);
    double *CB=mxGetPr(prhs[1]);
    int N=mxGetScalar(prhs[2]);
    int n=mxGetScalar(prhs[3]);

    float *llr=(float*)malloc(sizeof(float)*N);
    int *CBR=(int*)malloc(sizeof(int)*N);

    for(int i=0;i
        llr[i]=(float)a[i];
        CBR[i]=(int)CB[i];
    }

    int *u_d;
    u_d=decode(llr,CBR,N,n);
    plhs[0]=mxCreateDoubleMatrix(1,N, mxREAL);
    double *u=mxGetPr(plhs[0]);
    for(int i=0;i
        u[i]=u_d[i];
    }
}

matlab测试函数

clc
clear
addpath('GA/')
N=8;
K=N/2;
n=log2(N);
SNR=4;
sigma=10^(-SNR/20);
[channels, ~] = GA(sigma, N);
[~, channel_ordered] = sort(channels, 'descend');%降序
info_bits = sort(channel_ordered(1 : K), 'ascend');
frozen_bits = ones(N , 1);
frozen_bits(info_bits) = 0;%把1当成冻结比特
info_bits_logical = logical(mod(frozen_bits + 1, 2))';
CBR=double(info_bits_logical);
xx=CBR;
%随机产生的信源
x=randi([0,1],1,K);
for i=1:K
    xx(info_bits(i))=x(i);
end
xx
%编码
u=mod(xx*G(n),2);
%调制
u1=1-2*u;
%加噪
nn=randn(1,N);
nn1=nn*sigma;
%叠加
y=u1+nn1;
%转换为对数域
llr=2*y/(sigma^2);
%译码
u_d=demo(llr,CBR,N,n)
count=0;
for i=1:N
    if u_d(i)~=xx(i)
        count=count+1;
    end
end