天线选择定点算法

天线选择定点算法
一、 系统平台与原理算法间的对应关系
在文档[1]中提出，用以下两个公式来对不同天线配置下的信道容量进行遍历式计算，然后通过
比较这些容量的大小来确定天线选择方案：
()()1
2
,0
(_)log det N i M
i C m im o index M -===
∑A
（1）
2
H s i M M i i
E M σ
+
A =I H H ，%% （2） 其中i M A ，为Hermitian 矩阵，维数为M M ⨯（[]1,4M ∈）。

i
H %是第i 个参考信号子载波的有效信道矩阵，维数为4M ⨯。

s E 是总的信号的发送能量（归一化能量发送时取值为1），M 是发送流数，2
σ是多个接收天线上的平均噪声方差。

N 为参与运算的总的子载波信道数目。

在系统算法平台中，发射端的调制信号能量归一化以后，经过IFFT 运算后总能量衰减为原来的1/1024，总的信号发送能量没有进行归一化，在对过了多径的信号噪声的时候，噪声功率计算为：
_H POWER ： 多径信道增益 _E symbol ： 平均每个符号的能量，星座点归一化为1 T N ：
发射天线数
_FFT NOTES : FFT 点数
2
σ：
噪声功率谱密度 0
b E N ：
指定的信噪比
M od ：
调制方式
实际加入的噪声是信号的实部和虚部分别加入的，功率谱密度为：
2
__2
b H POW ER E symbol
E M od N σ
⨯=
⨯
⨯
在利用CP 或是preamble 进行2
σ的估计时，由于是在时域进行的，接收到的信号经过AGC 以后会有个AGC 的调整因子：
_agc divisor ：AGC 调整因子
这样估计到的噪声功率谱密度是2σ的倍数，倍数关系为：2_agc divisor ，即估计得到的噪声功率谱密度：
222_est agc divisor σσ=⨯
由于接收到的midamble 也是经过AGC 的，且信道估计是在频域进行，而接收机这端FFT 运算也没有归一化，会对信号有个放大，从仿真结果看放大倍数在_FFT NOTES 左右，没有作大量的统计仿真，观察了几组数据，选择用_FFT NOTES 来作为估算结果。

综合考虑发送端、接收端的IFFT 、FFT 和AGC 的影响，估计得到的频域信道和实际频域信道间的关系为：
i =H %%
以及算法中估算的结果，公式（2）变为以下形式：
2
2
2
2
___))_H s i M M i i
T
H M est H T M i i
E M E sym bol N FFT N O TES M E sym bol N M FFT N O TES
σ
σσ+
⨯=+⨯=+
⨯A =I H H I I (H (H ，%%%%%%（3）
二、
定点算法中各个模块及各变量位宽设置
有了上述结果，下面给出天线选择各个模块的定点化程序。

2.1 噪声功率谱密度计算
通过估计得到的噪声功率谱密度及其位宽设置为： 2
est σ 位宽为（18, 18），当
b E N 的变换范围为：0~30dB 的时候，估算得到的2
est σ的浮点数
值范围为（0.0000488281, 0.0488281）
通过查表法求得其倒数为：
21
est
σ
位宽为（18, 2）
2.2 天线选择模块
void P_rank_cal( P_Complex* h_ch, //输入的信道矩阵
int input_w,int input_w_d, //信道矩阵位宽（18,10）
int N, //天线选择中所选用的子载波数
FIX var_inverse, //信道矩阵位宽（18,10）
int var_w, int var_w_d, //噪声功率谱密度的倒数位宽（18,2）
int *mimo_indicator, //计算得到的mimo索引值
int nr) //接收天线数
{
int eff_nt,i,j,k,nr_idx,nt_idx,mimo_idx;
int h_eff_w=input_w; //计算有效信道位宽和输入总位宽一致，
int h_eff_w_d=input_w_d-2; //有效信道小数位宽减两位，分两位计入整数位宽
int matrix_a_w=18; //H’*H得到的矩阵的总位宽，统一为18
int matrix_a_w_d=4; // H’*H得到的矩阵的小数位宽，统一为4
int matrix_adjust2_a_w=18; //公式（2）中计算后的总位宽
int matrix_adjust2_a_w_d=6; //公式（2）中计算后的总小数位宽
int det_w=18; //算得的行列式位宽
int det_w_d=8; //算得的行列式小数位宽
FIX det=0;
FIX temp=0;
FIX tmp1,tmp2;
FIX* temp_capa_fix=(FIX*)malloc(sizeof(FIX)); //对每个载波计算算得的容量，最后要累加到变量capacity中，共32个累加量（N）int temp_capa_fix_w=11; //每个累加量的总位宽
int temp_capa_fix_w_d=7; //每个累加量的小数位宽
FIX capacity[4]; //现在只计算天线独立数为1和4时的容量，每个位宽为（16,7），由32个temp_capa_fix累加得到
int capacity_w=16;
int capacity_w_d=7;
/*****************************************************************************
所需变量空间分配，可以不管
*****************************************************************************/
P_Complex ***h_matrix,***h_eff,*h_subcarr,*h_subcarr_Htranspose, *matrixA, tmp;
h_eff = alloc_mem_P_complex3(N,nr,4);
h_subcarr = alloc_mem_P_complex1(nr*4);
h_subcarr_Htranspose = alloc_mem_P_complex1(4*nr);
matrixA = alloc_mem_P_complex1(4*4);
memset(capacity,0,sizeof(FIX)*4); // mimo_index < 4
for (mimo_idx=0;mimo_idx<4;mimo_idx+=3) // mimo_index < 4
{
/*利用函数P_eff_ch_cal和输入的h_ch来计算对应于不同mimo_idx下的有效信道*/
P_eff_ch_cal(h_ch, input_w, input_w_d,N, mimo_idx,h_eff,h_eff_w,h_eff_w_d, &eff_nt,nr);
for (k=0;k<N;k++)
{
/*每次提取一个载波上的有效信道矩阵到变量h_subcarr中,即公式（2）中的Hi*/
h_subcarr h_eff[k*nr_idx*nt_idx];
/*求Hi共轭转置 Hi’*/
P_Htranspose(h_subcarr,nr,eff_nt,h_subcarr_Htranspose,1);
/*矩阵相乘:Hi*Hi’*/
P_Matrix_Multiply(matrixA, h_subcarr_Htranspose,h_subcarr,h_eff_w,h_eff_w_d,matrix_a_w,matrix_a_w_d, eff_nt, eff_nt, nr);
/* 求行列式对应的矩阵,即公式（2）中的矩阵，由于有效天线数为1的时候，要对计算得到的矩阵除以2^18，而有效天线数为4时，要除以2^20，所以此处将分根有效天线数的不同，对（2）中矩阵的计算分开来处理，但是最终结果都是保存为（18，6）*/
if(eff_nt==1)
{
for (i=0;i<eff_nt;i++)
{
for (j=0;j<eff_nt;j++)
{
tmp1=Multiply((matrixA+i*eff_nt+j)->re,var_inverse);
tmp2=Multiply((matrixA+i*eff_nt+j)->im,var_inverse);
(matrixA+i*eff_nt+j)->re= AdjustFIX(tmp1>>18,36,6,matrix_adjust2_a_w,matrix_adjust2_a_w_d);
(matrixA+i*eff_nt+j)->im= AdjustFIX(tmp2>>18,36,6,matrix_adjust2_a_w,matrix_adjust2_a_w_d);
if (i==j)
(matrixA+i*eff_nt+j)->re = (matrixA+i*eff_nt+j)->re + Tofixed(1,matrix_adjust2_a_w,matrix_adjust2_a_w_d);
}
}
}
else if(eff_nt==4)
{
for (i=0;i<eff_nt;i++)
{
for (j=0;j<eff_nt;j++)
{
tmp1=Multiply((matrixA+i*eff_nt+j)->re,var_inverse);
tmp2=Multiply((matrixA+i*eff_nt+j)->im,var_inverse);
(matrixA+i*eff_nt+j)->re= AdjustFIX(tmp1>>20,36,6,matrix_adjust2_a_w,matrix_adjust2_a_w_d);
(matrixA+i*eff_nt+j)->im= AdjustFIX(tmp2>>20,36,6,matrix_adjust2_a_w,matrix_adjust2_a_w_d);
if (i==j)
(matrixA+i*eff_nt+j)->re = (matrixA+i*eff_nt+j)->re + Tofixed(1,matrix_adjust2_a_w,matrix_adjust2_a_w_d);
}
}
}
/* 求行列式*/
det = P_Matrix_Determinant(matrixA, eff_nt,matrix_adjust2_a_w,matrix_adjust2_a_w_d,det_w,det_w_d);
if(det.re < Tofixed(1,det_w,det_w_d)) // calculation protection
{
det.re = Tofixed(1,det_w,det_w_d);
}
*temp_capa_fix=0;
/* 查表法求对数*/
LOG2FORFIX(tmp.re, det_w,det_w_d, temp_capa_fix, temp_capa_fix_w, temp_capa_fix_w_d);
/* 计算得到的容量累积*/
capacity[mimo_idx]= Add(AdjustFIX(*temp_capa_fix,temp_capa_fix_w, temp_capa_fix_w_d,capacity_w,capacity_w_d),capacity[mimo_idx]);
}
}
/* 选择最大容量的天线配置为系统下一帧的天线数*/
for (mimo_idx=0;mimo_idx<4;mimo_idx++)
{
if (capacity[mimo_idx]>temp)
{
temp = capacity[mimo_idx];
*mimo_indicator= mimo_idx;
}
}
}
2.3 各个子模块
2.3.1 有效信道计算模块
根据输入的mimo_index值来计算有效信道
输入、输出参数位宽参见上一节。

不管信道是几维的，4*4还是1*1，最终都将位宽调整为（18,8）。

void P_eff_ch_cal(P_Complex ***h_ch, int input_w,int input_w_d,int N, int mimo_index, P_Complex ***h_eff,int h_eff_w,int h_eff_w_d, int *eff_nt,int nr)
{
int nt_idx,nr_idx,k;
P_Complex t1,t2;
/*分别根据输入的mimo_index不同，对频域信道矩阵做不同处理*/
if (mimo_index==0)
{
for (nt_idx=0;nt_idx<4;nt_idx++)
{
for (nr_idx=0;nr_idx<nr;nr_idx++)
{
for (k=0;k<N;k++)
{
h_eff[k][nr_idx][nt_idx] = AdjustPC(h_ch[k][nr_idx][nt_idx],input_w, input_w_d,h_eff_w, h_eff_w_d);
}
}
}
*eff_nt=4;
}
//mimo_index为1时，0、1两根天线上的数据相加
else if (mimo_index==1)
{
for (nr_idx=0;nr_idx<nr;nr_idx++)
{
for (k=0;k<N;k++)
{ h_eff[k][nr_idx][0] = P_Complex_Add( h_ch[k][nr_idx][0] , h_ch[k][nr_idx][1] );
h_eff[k][nr_idx][0] = AdjustPC(h_eff[k][nr_idx][0],input_w+1,input_w_d,h_eff_w, h_eff_w_d);
h_eff[k][nr_idx][1] = AdjustPC(h_ch[k][nr_idx][2],input_w,input_w_d,h_eff_w, h_eff_w_d);
h_eff[k][nr_idx][2] = AdjustPC(h_ch[k][nr_idx][3],input_w,input_w_d,h_eff_w, h_eff_w_d);
}
}
*eff_nt=3;
}
//mimo_index为1时，0,1两根天线上的数据相加，2、3两根天线上的数据相加
else if (mimo_index==2) // mimo_index = 7
{
for (nr_idx=0;nr_idx<nr;nr_idx++)
{
for (k=0;k<N;k++)
{
h_eff[k][nr_idx][0] = P_Complex_Add( h_ch[k][nr_idx][0] , h_ch[k][nr_idx][1] );
h_eff[k][nr_idx][0] = AdjustPC(h_eff[k][nr_idx][0],input_w+1,input_w_d,h_eff_w, h_eff_w_d);
h_eff[k][nr_idx][1] = P_Complex_Add( h_ch[k][nr_idx][2] , h_ch[k][nr_idx][3] );
h_eff[k][nr_idx][1] = AdjustPC(h_eff[k][nr_idx][1],input_w+1,input_w_d,h_eff_w, h_eff_w_d); }
}
*eff_nt=2;
}
//mimo_index为1时，0,1， 2、3两根天线上的数据相加
else if (mimo_index==3) // mimo_index = 10
{
for (nr_idx=0;nr_idx<nr;nr_idx++)
{
for (k=0;k<N;k++)
{
t1 = P_Complex_Add( h_ch[k][nr_idx][0] , h_ch[k][nr_idx][1] );
t2 = P_Complex_Add( h_ch[k][nr_idx][2] , h_ch[k][nr_idx][3] );
h_eff[k][nr_idx][0] = P_Complex_Add(t1, t2 );
h_eff[k][nr_idx][0] = AdjustPC(h_eff[k][nr_idx][0],input_w+2,input_w_d,h_eff_w, h_eff_w_d); }
}
*eff_nt=1;
}
else
{
printf("Wrong mimo_index in eff_ch_cal.c\n");
}
return;
}
matrixA中
void P_Matrix_Multiply(P_Complex * matrixA,
P_Complex * matrixB,
P_Complex * matrixC,
int input_w,
int input_w_d,
int matrix_a_w,
int matrix_a_w_d,
int row,
int column,
int same_parameter
)
{
//FILE* fptest=fopen("f:\\fix.txt","at+");
int i,j,k;
P_Complex * move_pointerA;
const P_Complex * move_pointerB;
const P_Complex * move_pointerC;
memset(matrixA,0,sizeof(P_Complex)*row*column);
for(i=0;i<row;i++)
{
move_pointerB=matrixB+i*same_parameter;
for(j=0;j<column;j++)
{
move_pointerA=matrixA+i*column+j;
move_pointerC=matrixC+j;
for(k=0;k<same_parameter;k++)
{
*move_pointerA=P_Complex_Add(*move_pointerA,P_Complex_Multiply(*(move_pointerB+k),*(move_pointerC+k*column)));
}
*move_pointerA=AdjustPC(*move_pointerA,input_w*2+2,input_w_d*2,matrix_a_w,matrix_a_w_d);
}
}
}
行列式是用递归的方法来求，递归过程中所有中间结果位宽统一调整为（18,8）
P_Complex P_Matrix_Determinant(P_Complex * matrix,
int eff_nt,
int input_w,
int input_w_d,
int det_w,
int det_w_d)
{
P_Complex * matrixA = matrix;
P_Complex determinant;
//根据有效信道数来确定矩阵的维数
switch (eff_nt)
{
case 1:
determinant = AdjustPC(*matrixA,input_w,input_w_d,det_w,det_w_d);
break;
case 2:
determinant = P_Matrix_Determinant_Two(matrixA,input_w,input_w_d,det_w,det_w_d);
break;
case 3:
determinant = P_Matrix_Determinant_Three(matrixA, input_w,input_w_d,det_w,det_w_d);
break;
case 4:
determinant = P_Matrix_Determinant_Four(matrixA, input_w,input_w_d,det_w,det_w_d);
}
return determinant;
}
// N = column 二阶行列式,对角线相乘得到
P_Complex P_Matrix_Determinant_Two(P_Complex * matrix,
int input_w,
int input_w_d,
int det_w,
int det_w_d)
{
P_Complex determinant = {0,0};
P_Complex temp1,temp2;
temp1 = P_Complex_Multiply(*(matrix), *(matrix+3));
temp2 = P_Complex_Multiply(*(matrix+1), *(matrix+2));
determinant = P_Complex_Minus( temp1 , temp2 );
// 调整行列式位宽为啊（18,8）
determinant = AdjustPC(determinant,input_w*2,input_w_d*2,det_w,det_w_d);
return(determinant);
}
// N = column 三阶行列式,利用二阶行列数递归得到
P_Complex P_Matrix_Determinant_Three(P_Complex * matrix,
int input_w,
int input_w_d,
int det_w,
int det_w_d)
{
int i;
int N=3;
P_Complex tmp;
P_Complex determinant = {0,0};
P_Complex * matrix_part;
matrix_part = (P_Complex*)malloc((N-1)*(N-1)*sizeof(P_Complex));
for(i=0; i<N; i++)
{
// 先求代数余子式
P_Matrix_Part(matrix_part, matrix, N, N, 0, i);
P_Scalar_Multiply2(matrix_part,input_w,input_w_d, N-1, N-1, i);
// 求得的代数余子式和元素相乘
tmp= P_Complex_Multiply(*(matrix+i), P_Matrix_Determinant_Two(matrix_part,input_w,input_w_d,det_w,det_w_d));
tmp=AdjustPC(tmp,input_w+det_w,input_w_d+det_w_d,det_w,det_w_d);
determinant = P_Complex_Add( determinant, tmp);
}
// 调整行列式位宽为啊（18,8）
determinant = AdjustPC(determinant,det_w+2,det_w_d,det_w,det_w_d);
free(matrix_part);
return(determinant);
}
// N = column 四阶行列式,利用三阶行列数递归得到
P_Complex P_Matrix_Determinant_Four(P_Complex * matrix,
int input_w,
int input_w_d,
int det_w,
int det_w_d)
{
int i;
int N=4;
P_Complex tmp;
P_Complex determinant = {0,0};
P_Complex *matrix_part;
matrix_part = (P_Complex *) malloc((N-1)*(N-1)*sizeof(P_Complex));
for(i=0; i<N; i++)
{
// 先求代数余子式
P_Matrix_Part(matrix_part, matrix, N, N, 0, i);
P_Scalar_Multiply2(matrix_part,input_w,input_w_d, N-1, N-1, i);
// 求得的代数余子式和元素相乘
tmp= P_Complex_Multiply(*(matrix+i), P_Matrix_Determinant_Three(matrix_part, input_w,input_w_d,det_w,det_w_d));
tmp=AdjustPC(tmp,input_w+det_w,input_w_d+det_w_d,det_w,det_w_d);
determinant = P_Complex_Add( determinant, tmp);
}
// 调整位宽
determinant = AdjustPC(determinant,det_w+2,det_w_d,det_w,det_w_d);
free(matrix_part);
return(determinant);
}
2.3.4 查表法求对数
算法由朱正航提供
利用查表法求对数，表长为2^7+1=129,表可事先生成，生成以后保存起来，其中表中的数据位宽为（8,7），即小数位7位，生成函数如下：
void LOG_TABLE_generate(FIX* log_table, int table_length)
{
int i=0;
for(i=0;i<table_length;i++)
{
log_table[i]=Tofixed(log10(1+i*0.0078125)/log10(2),8,7);
}
}
利用上面生成的表格，对输入的数据data_in进行查表，得到对数值。

找到第一个“1”以后，其后7位数据有效，以该7位数据为地址直接索引表格对应的定点对数值，再和第一个“1”相对于小数点的位置值相加得到对数值。

void LOG2FORFIX(FIX data_in, int input_w, int input_w_d, FIX* log_value, int log_value_w, int log_value_w_d)
{
int i=0,j=0;
char table_index=0;
char location=input_w-input_w_d-1;
FIX* log_table=(FIX*)malloc(sizeof(FIX)*129); //3bits
FIX tmp=Tofixed(256,input_w,input_w_d);
// 生成对数表格，
LOG_TABLE_generate(log_table,129);
// 寻找第一个“1”所在的位置：location（相对于小数点），并找到其后7位有效数据，存入table_index，用于查表
for(i=0;i<input_w;i++)
{
// 找到“1”以后，计算table_index
if((tmp & data_in)!=0)
{
for(j=0;j<7;j++)
{
tmp=tmp>>1;
table_index=table_index << 1;
if((tmp & data_in) !=0)
{
table_index=table_index | 1;
}
}
// 确定第一个“1”相对于小数点的位置
location=location-i-1;
break;
}
else
tmp=tmp>>1;
}
// location+table[table_index]得到最后的对数值
*log_value=AdjustFIX(log_table[table_index],8,7,log_value_w,log_value_w_d)+Tofixed((float)location,log_value_w, log_value_w_d);
free(log_table);
}
// 矩阵数乘：确定代数余子式，系数（-1）^(i+j)和子矩阵第一行相乘
void P_Scalar_Multiply2(P_Complex * matrixA,
int input_w,
int input_w_d,
int row,
int column,
int scalar
)
{
int i,j;
// 由于乘数为±1，所以只需要根据输入的scalar来确定是不变还是取反即可
if(scalar%2==0)
{
for(i=0; i<1; i++)
{
for(j=0; j<column; j++)
{
(matrixA + i*column + j)->re = (matrixA + i*column + j)->re;
(matrixA + i*column + j)->im = (matrixA + i*column + j)->im;
}
}
}
else
{
for(i=0; i<1; i++)
{
for(j=0; j<column; j++)
{
(matrixA + i*column + j)->re = -(matrixA + i*column + j)->re;
(matrixA + i*column + j)->im = -(matrixA + i*column + j)->im;
}
}
}
}
// 除去第i行，第j列后剩下的矩阵
void P_Matrix_Part(P_Complex * matrixA,
P_Complex * matrixB,
int row,
int column,
int I,
int J)
{
int i,j;
P_Complex * move_pointerA ;
move_pointerA = matrixA;
for(i=0; i<row; i++)
{
for(j=0; j<column; j++)
{
if(i == I || j == J)
continue;
else
{
move_pointerA->re = (matrixB + i*column + j)->re;
move_pointerA->im = (matrixB + i*column + j)->im;
move_pointerA++;
}
}
}
}。