📄 ofdm.c
字号:
// ---------------------------------------------
// 参数定义
// ---------------------------------------------
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <complex>
// 变量初始化
int NumLoop = 500;
int NumSubc = 128;
int NumCP = 8;
int SyncDelay = 0;
//-----------------------------------------------
// 子载波数 128
// 位数/ 符号 2
// 符号数/ 载波 500
// 训练符号数 0
// 循环前缀长度 8 (1/16)*T
// 调制方式 4-QAM
// 多径信道数 3
// IFFT Size 128
// 信道最大时延 2
//-----------------------------------------------
// ---------------------------------------------
// QAM 调制
// ---------------------------------------------
// 产生随机二进制输入数据
void
data_in()
{
long int p;
int i;
char datatmp=0;
int BitsTx[64000]={0};
for (i=0;i<NumLoop*NumSubc;i++)
{
p=rand();
if (p/16383<1) {datatmp=0;} //p是随机数，如果小于阈值（式中为16383）
else
{
datatmp=1;
}
BitsTx[i]=datatmp;
}
}
//将二进制数据转化为极性数据
//void
//bin2pol()
//{
// for (i=0;i<64000;i++)
// if (datatmp[i]==0)
// data[i]=-1;
// else
// data[i]=1;
//}
//调整数组
void
data_reshape()
{
int BitsTx[64000];
int QAM_tmp[2][32000]={0};
int i;
for(j=0;j<32000;j++)
{
for(k=0;k<2;k++)
{
QAM_tmp[k][j]=BitsTx[i];
i++;
}
}
// 调制 (产生 QAM 符号->s/p)
// 输入: BitsTx[NumLoop*NumSubc]; 输出: SymQAM(NumLoop,NumSubc/2)
//void
//data_QAM()
//{
// int BitsTxtmp1[32000]={0};
// int BitsTxtmp2[32000]={0};
// int BitsTxtmp[128][500];
// int j,k;
//int i=0;
// int BitsTx[64000];
// for (k=0;k<128;k++)
// {
// for (j=0;j<500;j++)
// BitsTxtmp[k][j]=BitsTx[i];
// i++;
// }
//}
//SymQAMtmptmp = bi2de(SymQAMtmp,2,'left-msb');?
// QAM 调制
void
QAM_mod()
{
int k;
int j=0;
int Bits_tmp[2][32000];
typedef complex<float> SymQAM[32000];
for (k=0;k<32000;k++)
{
if
((BitsTxtmp[j][k]==0)&(BitsTxtmp[j+1][k]==0))
SymQAM[k]=1+i;
if
((BitsTxtmp[j][k]==0)&(BitsTxtmp[j+1][k]==1))
SymQAM[k]=-1+i;
if
((BitsTxtmp[j][k]==1)&(BitsTxtmp[j+1][k]==0))
SymQAM[k]=-1-i;
if
((BitsTxtmp[j][k]==1)&(BitsTxtmp[j+1][k]==1))
SymQAM[k]=1-i;
}
}
// ---------------------------------------------
// IFFT
// ---------------------------------------------
// 输入: SymQAM[NumLoop][NumSubc/2]; 输出: SymIFFT[NumSubc][NumLoop]
int l=1;
int n=32000;
void
Sym_fft(double ifft_in_real[], double ifft__in_imag[], int n, int k, double ifft_out_real[], double ifft_out_image[], int l)
{
// 入口参数：
// l: l = 0, 傅立叶变换; l = 1, 逆傅立叶变换
// il: il = 0,不计算傅立叶变换或逆变换模和幅角；il = 1,计算模和幅角
// n: 输入的点数，为偶数，一般为32，64，128，...,1024等
// k: 满足n=2^k(k>0),实质上k是n个采样数据可以分解为偶次幂和奇次幂的次数
// ifft_in_real[]: l=0时，存放N点采样数据的实部
// l=1时, 存放傅立叶变换的N个实部
// double ifft__in_imag[]: l=0时，存放N点采样数据的虚部
// l=1时, 存放傅立叶变换的N个虚部
//
// 出口参数：
// ifft_out_real: l=0, 返回傅立叶变换的实部
// l=1, 返回逆傅立叶变换的实部
// ifft_out_image: l=0, 返回傅立叶变换的虚部
// l=1, 返回逆傅立叶变换的虚部
// pr[]: il = 1,i = 0 时，返回傅立叶变换的模
// il = 1,i = 1 时，返回逆傅立叶变换的模
// pi[]: il = 1,i = 0 时，返回傅立叶变换的辐角
// il = 1,i = 1 时，返回逆傅立叶变换的辐角
int it,m,is,i,j,nv,l0;
double p,q,s,vr,vi,poddr,poddi;
for (it=0;it<=n-1;it++)
{
m=it;
is=0;
for(i=0;i<=k-1;i++)
{
j=m/2;
is=2*is+(m-2*j);
m=j;
}
fr[it]=pr[is];
fi[it]=pi[is];
}
//----------------------------
ifft_in_real[0]=1.0;
ifft__in_imag[0]=0.0;
p=6.283185306/(1.0*n);
pr[1]=cos(p);
pi[1]=-sin(p);
if (l!=0)
pi[1]=-pi[1];
for (i=2;i<=n-1;i++)
{
p=pr[i-1]*pr[1];
q=pi[i-1]*pi[1];
s=(pr[i-1]+pi[i-1])*(pr[1]+pi[1]);
pr[i]=p-q;
pi[i]=s-p-q;
}
for (it=0;it<=n-2;it=it+2)
{
vr=fr[it];
vi=fi[it];
fr[it]=vr+fr[it+1];
fi[it]=vi+fi[it+1];
fr[it+1]=vr-fr[it+1];
fi[it+1]=vi-fi[it+1];
}
m=n/2;
nv=2;
for (l0=k-2;l0>=0;l0--)
{
m=m/2;
nv=2*nv;
for(it=0;it<=(m-1)*nv;it=it+nv)
for (j=0;j<=(nv/2)-1;j++)
{
p=pr[m*j]*fr[it+j+nv/2];
q=pi[m*j]*fi[it+j+nv/2];
s=pr[m*j]+pi[m*j];
s=s*(fr[it+j+nv/2]+fi[it+j+nv/2]);
poddr=p-q;
poddi=s-p-q;
fr[it+j+nv/2]=fr[it+j]-poddr;
fi[it+j+nv/2]=fi[it+j]-poddi;
fr[it+j]=fr[it+j]+poddr;
fi[it+j]=fi[it+j]+poddi;
}
}
if(l!=0)
for(i=0;i<=n-1;i++)
{
fr[i]=fr[i]/(1.0*n);
fi[i]=fi[i]/(1.0*n);
}
if(il!=0)
for(i=0;i<=n-1;i++)
{
pr[i]=sqrt(fr[i]*fr[i]+fi[i]*fi[i]);
if(fabs(fr[i])<0.000001*fabs(fi[i]))
{
if ((fi[i]*fr[i])>0)
pi[i] = 90.0;
else
pi[i] = -90.0;
}
else
pi[i] = atan(fi[i]/fr[i])*360.0/6.283185306;
}
return;
}
//char SymIFFT[NumSubc][NumLoop] = 0;
//SymIFFTtmp = reshape(SymQAM,NumSubc/2,NumLoop);
//SymIFFTtmptmp = zeros(NumSubc,NumLoop);
//SymIFFTtmptmp(1,:) = real(SymIFFTtmp(1,:)); % 实数
//SymIFFTtmptmp(NumSubc/2+1,:) = imag(SymIFFTtmp(1,:)); % 实数
// 这么安排矩阵的目的是为了构造共轭对称矩阵
// 共轭对称矩阵的特点是 在ifft/fft的矢量上 N点的矢量
// 在0，N/2点必须是实数 一般选为0
// 1至N/2点 与 （N/2）+1至N-1点关于N/2共轭对称
//SymIFFTtmptmp(2:NumSubc/2,:) = SymIFFTtmp(2:NumSubc/2,:);
//SymIFFTtmptmp((NumSubc/2+2):NumSubc,:) = flipdim(conj(SymIFFTtmp(2:NumSubc/2,:)),1);
//%--------------------------------------------------------------------------
//% 函数说明：
//% B = flipdim(A,dim) returns A with dimension dim flipped.
//%/ When the value of dim is 1, the array is flipped row-wise down. When dim is 2,
//% the array is flipped columnwise left to right. flipdim(A,1) is the same as
//% flipud(A), and flipdim(A,2) is the same as fliplr(A).
//%--------------------------------------------------------------------------
//% >> a = [1 2 3; 4 5 6; 7 8 9; 10 11 12]
//% a =
//% 1 2 3
//% 4 5 6
//% 7 8 9
//% 10 11 12
//% >> b = flipdim(a,1)
//% b =
//% 10 11 12
//% 7 8 9
//% 4 5 6
//% 1 2 3
//SymIFFT = ifft(SymIFFTtmptmp,NumSubc,1);
// ---------------------------------------------
// 循环前缀
// ---------------------------------------------
// 输入: SymIFFT[NumSubc][NumLoop]; 输出: SymCP[NumAddPrefix][NumLoop]
void
Data_CP()
{
int NumAddPrefix;
int RowPrefix;
NumAddPrefix=NumSubc+NumCP;
int SymCP[NumAddPrefix][NumLoop]=0;
for (RowPrefix=NumSubc-NumCP+1;RowPrefix<=NumSubc;RowPrefix++;)
SymCP[NumAddPrefix][NumLoop]=SymIFFT[RowPrefix][NumLoop];
}
// ---------------------------------------------
// 通过信道
// ---------------------------------------------
// input: SymCP(NumSubc + NumCP,NumLoop); output: SymCh(1,(NumSubc + NumCP)*NumLoop)
// 输入：SymCP[NumAddPrefix][NumLoop]；输出：SymCP[1][NumAddPrefix*NumLoop]
//SymCh = zeros(1,(NumSubc + NumCP)*NumLoop);
void
Data_ch()
{
int SymCh[1][NumAddPrefix*NumLoop]=0;
int SymChtmp[1][NumAddPrefix*NumLoop]=0;
int SymChtmptmp[1][NumAddPrefix*NumLoop]=0;
int i;
int j=1;
for (i=1;i<=NumAddPrefix*NumLoop;i++)
{
SymChtmp[j][1]=SymCP[1][NumAddPrefix*NumLoop];
j++;
}
//SymChtmp = SymCP(:).'; % 进行这个转置操作之后就成了一个矢量
// % 相当于把矩阵的列向量依次排列 改变为一个行向量
//Ch = [1 1/2 1/4];
Ch=[1 1/2 1/4];
}
//SymChtmptmp = filter(Ch,1,SymChtmp);
//%--------------------------------------------------------------------------
//% 函数说明：
//% Firlter data with an infinite impulse response (IIR) or finite impulse response
//% (FIR) filter
//% y = filter(b,a,X) filters the data in vector X with the filter described by
//% numerator coefficient vector b and denominator coefficient vector a. If a(1) is
//% not equal to 1, filter normalizes the filter coefficients by a(1). If a(1) equals
//% 0, filter returns an error.
//%--------------------------------------------------------------------------
//% If X is a matrix, filter operates on the columns of X. If X is a multidimensional
//% array, filter operates on the first nonsingleton dimension.
//%-------------------------------------------------------------------------
// 添加AWGN?
//BerSnrTable = zeros(20,3);
//for snr=0:19; % = SNR + 10*log10(log2(2));
// BerSnrTable(snr+1,1) = snr;
//SymCh = awgn(SymChtmptmp,snr,'measured');
//%--------------------------------------------------------------------------
//% 函数说明：
//% AWGN Add white Gaussian noise to a signal.
//% Y = AWGN(X,SNR) adds white Gaussian noise to X. The SNR is in dB.
//% The power of X is assumed to be 0 dBW. If X is complex, then
//% AWGN adds complex noise.
//% ------------------------------------------------------------------------
//%
//% Y = AWGN(X,SNR,SIGPOWER) when SIGPOWER is numeric, it represents
//% the signal power in dBW. When SIGPOWER is 'measured', AWGN measures
//% the signal power before adding noise.
//%U.p5h.{.f8M.b0% -------------------------------------------------------------------------
//% Y = AWGN(X,SNR,SIGPOWER,STATE) resets the state of RANDN to STATE.
//%
//% Y = AWGN(..., POWERTYPE) specifies the units of SNR and SIGPOWER.
//%i.?#B0T0% POWERTYPE can be 'db' or 'linear'. If POWERTYPE is 'db', then SNR
//% is measured in dB and SIGPOWER is measured in dBW. If POWERTYPE is
//% 'linear', then SNR is measured as a ratio and SIGPOWER is measured
//% in Watts.
//%
//% Example: To specify the power of X to be 0 dBW and add noise to produce
//% an SNR of 10dB, use:
//% X = sqrt(2)*sin(0:pi/8:6*pi);
//% Y = AWGN(X,10,0);
//%
//% Example: To specify the power of X to be 0 dBW, set RANDN to the
//% 1234th
//% state and add noise to produce an SNR of 10dB, use:
//% X = sqrt(2)*sin(0:pi/8:6*pi);
//% Y = AWGN(X,10,0,1234);
//%
//% Example: To specify the power of X to be 3 Watts and add noise to
//% produce a linear SNR of 4, use:
//% X = sqrt(2)*sin(0:pi/8:6*pi);
//% Y = AWGN(X,4,3,'linear');
//% Example: To cause AWGN to measure the power of X, set RANDN to the
//% 1234th state and add noise to produce a linear SNR of 4, use:
//% X = sqrt(2)*sin(0:pi/8:6*pi);
//% Y = AWGN(X,4,'measured',1234,'linear');
// ---------------------------------------------
// 移除保护间隔
// ---------------------------------------------
// input: SymCh(1,(NumSubc + NumCP)*NumLoop); output: SymDeCP(NumSubc,NumLoop)
// 输入：SymCh[NumAddPrefix*NumLoop]；输出：SymDeCP[NumSubc][NumLoop]
//SymDeCP = zeros(NumSubc,NumLoop);
void
Data_deCP()
{
int SymDeCP[NumSubc][NumLoop]=0;
int SymDeCPtmp[NumSubc][NumLoop]=0;
int i,j;
int k=0;
for (i=1;i<=NumLoop;i++)
{
for(j=1;j<=NumAddPrefix;j++)
{
SymDeCPtmp[j][i]=SymCh[k];
k++;
}
SymDeCp[NumSubc][NumLoop]=SymDeCPtmp[NumSubc][Loop];
}
//SymDeCPtmp = reshape(SymCh,NumAddPrefix,NumLoop);
//SymDeCP = SymDeCPtmp((NumCP+1+SyncDelay):NumAddPrefix+SyncDelay,:);
// ---------------------------------------------
// FFT
// ---------------------------------------------
/* 采样来的数据放在dataR[ ]数组中，运算前dataI[ ]数组初始化为0 */
void FFT(float dataR[],float dataI[])
{
int x0,x1,x2,x3,x4,x5,x6;
int L,j,k,b,p;
float TR,TI,temp;
/********** following code invert sequence ************/
for(i=0;i<128;i++)
{
x0=x1=x2=x3=x4=x5=x6=0;
x0=i&0x01;
x1=(i/2)&0x01;
x2=(i/4)&0x01;
x3=(i/8)&0x01;
x4=(i/16)&0x01;
x5=(i/32)&0x01;
x6=(i/64)&0x01;
xx=x0*64+x1*32+x2*16+x3*8+x4*4+x5*2+x6;
dataI[xx]=dataR[i];
}
for(i=0;i<128;i++)
{
dataR[i]=dataI[i];
dataI[i]=0;
}
/************** following code FFT *******************/
for(L=1;L<=7;L++)
{/* for(1) */
b=1;
i=L-1;
while(i>0)
{
b=b*2;
i--;
}/*b=2^(L-1)*/
for(j=0;j<=b-1;j++)/*for(2)*/
{
p=1;
i=7-L;
while(i>0)/*p=pow(2,7-L)*j;*/
{
p=p*2;
i--;
}
p=p*j;
for(k=j;k<128;k=k+2*b)/*for(3)*/
{
TR=dataR[k];
TI=dataI[k];
temp=dataR[k+b];
dataR[k]=dataR[k]+dataR[k+b]*cos_tab[p]+dataI[k+b]*sin_tab[p];
dataI[k]=dataI[k]-dataR[k+b]*sin_tab[p]+dataI[k+b]*cos_tab[p];
dataR[k+b]=TR-dataR[k+b]*cos_tab[p]-dataI[k+b]*sin_tab[p];
dataI[k+b]=TI+temp*sin_tab[p]-dataI[k+b]*cos_tab[p];
} /* END for (3) */
} /* END for (2) */
} /* END for (1) */
for(i=0;i<32;i++){ /* 只需要32次以下的谐波进行分析 */
w[i]=sqrt(dataR[i]*dataR[i]+dataI[i]*dataI[i]);
w[i]=w[i]/64
}
w[0]=w[0]/2;
} /* END FFT */
//SymFFT = fft(SymDeCP,NumSubc,1);
// ---------------------------------------------
// 判决(包括QAM解码)
// ---------------------------------------------
// SymFFT(NumSubc,NumLoop); output: SymDec(NumSubc,NumLoop)
// 输入：SymFFT[NumSubc][NumLoop];输出：SymDec[NumSubc][NumLoop]
//SymDec = zeros(NumSubc,NumLoop);
void
Data_Dec()
{
int SymDec[NumSubc][NumLoop]=0;
//SymEqtmp(1,:) = SymFFT(1,:)+i*SymFFT(NumSubc/2+1,:);?
//SymEqtmp(2:NumSubc/2,:) = SymFFT(2:NumSubc/2,:);?
for (m=1;m<=NumLoop;m++)
{
for (n=1;n<=NumSubc/2;n++)
{
Real = real(SymEqtmp(n,m));
Imag = imag(SymEqtmp(n,m));
if( abs((Real -1)) < abs((Real +1)))
SymDec[2*n-1][m]=1;
else
SymDec[2*n-1][m]=0;
if( abs((Imag -1)) < abs((Imag +1 )) )
SymDec[2*n][m]=1;
else
SymDec[2*n][m]=0;
}
}
// 误码率
void
Data_err()
{
int BitsRx[NumSubc*NumLoop]=0;
int i;
int j=1;
int k=1;
int m=0;
int n;
for (i=1;i<=NumLoop;i++)
{
for(k=1;k<=NumSubc;k++)
{
BitsRx[j] = SymDec[k][i];
j++;
}
for (j=1;j<=NumSubc*NumLoop;j++)
if (BitsTx[j]!=BitsRx)
m++;
n=m/(NumSubc*NumLoop)*100%;
printf('%d',m);
printf('%d',n);
}
// ---------------------------------------------
// The END
// ---------------------------------------------
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -