%OFDM程序示例
%信道选用 J.G.Proakis Digital Communications 影印版 第四版 P631
%采用64子信道,应用16QAM传输,每个子载波均传输相同的比特数4。
%不使用信道编码;使用迫零均衡和MMSE均衡,比较两者性能。
%试应用仿真的方法画出性能曲线。
close all;
clear all;
N=64; %The number of carriers
%h=[1];
%h=[0.407,0.815,0.407];
h=[0.04,-0.05,0.07,-0.21,-0.5,0.72,0.36,0,0.21,0.03,0.07]; %channel impulse response
L=length(h);
v=length(h)-1;%the size of the cyclic prefix
b=4; %For 16QAM, one symbol carries source bits
C16QAM=[3+3j,1+3j,-1+3j,-3+3j,-3+j,-1+j,1+j,3+j,3-j,1-j,-1-j,-3-j,-3-3j,-1-3j,1-3j,3-3j];
B16QAM=[0 0 0 0;0 0 0 1;0 0 1 1;0 0 1 0;0 1 1 0;0 1 1 1;0 1 0 1;0 1 0 0;
1 1 0 0;1 1 0 1;1 1 1 1;1 1 1 0;1 0 1 0;1 0 1 1;1 0 0 1;1 0 0 0];
srcBlockSize=N*b;
q=2^b; %the size of the signal set
Es=real(sum(C16QAM.*conj(C16QAM))/q);%Es= sum((abs(C16QAM)).^2)/q;
Es=Es/N; % If noise is generated in time domain, Es should be multiplied by a factor 1/N
Eb=Es/b;
Na=N+v;%The size of the transmitted block with cyclic prefix
Nb=N+v+L-1;%The size of the liner convolutional result of the transmitted block and the channel
noiseVec=zeros(1,Nb);
X=zeros(1,N); % one source symbol block
x=zeros(1,N); % for the IFFT of X
xt=zeros(1,Na); % transmitted time domain block with cyclic prefix
yt=zeros(1,Nb);%received time domain block
y=zeros(1,N);%received time domain signal removal of CP
Y=zeros(1,N);%the vector fot FFT of y
XD=zeros(1,N);
YD=zeros(1,N);
dist=zeros(1,q);
dist1=zeros(1,q);
src=zeros(1,srcBlockSize);
rev=zeros(1,srcBlockSize);
index=0;
index1=0;
H=fft(h,N); %the frequency domain character of the channel
randn('state',sum(100*clock)); %Initialize RANDN to a different state each time.
recIndex=1;%an index for recoder buffer
recIndex1=1;
tic
for EbN0=0:1:20
N0=Eb*10^(-EbN0/10);
noiseRoot=sqrt(N0/2);
errorCount=0;
errorCount1=0;
testLength=0;
BER=0;
BER1=0;
while(1)
testLength=testLength+srcBlockSize;
src=(randn(1,srcBlockSize)>0);% Generate a source block (equal probability)
%The following part complement the conversion of binary to q-QAM symbol
%a faster method is to use an index table,here we convert binary to symbol by calculating
for k=1:1:N
startIndex=(k-1)*b+1;
endIndex=startIndex+b-1;
vec=src(startIndex:endIndex);
for t=1:1:q
temp=sum(abs(vec-B16QAM(t,:)));
if(temp<1e-5)
index=t;
break;
end
end
X(k)=C16QAM(index);
end
%Mapping for frequency domain to time domain (multi-carrier modulation)
x=ifft(X); %IFFT
%Add CP
xt(1:v)=x( (N-v+1):N);
xt((v+1):(N+v))=x(1:N);
%Passing the channel with noise
yt=conv(xt,h);
noiseVec=randn(1,Nb)+j*randn(1,Nb);%Generate noise
noiseVec=noiseVec*noiseRoot;
yt=yt+noiseVec;
%Remove the cyclic prefix
y=yt((v+1):(v+N));
%FFT of y
Y=fft(y);
%Equalizer on the frequency domain
XD=Y./H; %Zero forcing equalizing ,alse Ls
YD= Y.*((conj(H))./((H.*conj(H))+mean(noiseVec.*conj(noiseVec))/mean(x.*conj(x))));%MMSE forcing equalizing
for k=1:1:N
%minminum distance decision
dist=abs(XD(k)-C16QAM);
dist1=abs(YD(k)-C16QAM);
[temp,index]=min(dist);
[temp,index1]=min(dist1);
vec=B16QAM(index,:);
vec1=B16QAM(index1,:);
startIndex=(k-1)*b+1;
endIndex=startIndex+b-1;
rev(startIndex:endIndex)=vec;
rev1(startIndex:endIndex)=vec1;
end
errorNum=sum(xor(src,rev));
errorNum1=sum(xor(src,rev1));
errorCount=errorCount+errorNum;
errorCount1=errorCount1+errorNum1;
BER=errorCount/testLength;
BER1=errorCount1/testLength;
% plot(real(XD),imag(XD),'*'); %This can be used to watch the received signal set
if((BER>0)&&(BER1>0))
if(EbN0<=9)
temp=200/BER;
temp1=200/BER1;
else
temp=100/BER;
temp1=100/BER1;
end
if((testLength>temp)&&(testLength>temp1))
if(testLength>150000)
break;
end
end
end
if(testLength>90000000)
break;
end
end
BER_Rec(recIndex)=BER
BER_Rec1(recIndex)=BER1
EbN0_Rec(recIndex)=EbN0
testLength_Rec(recIndex)=testLength;
recIndex=recIndex+1;
end
toc
%Save the result
save ofdm_data EbN0_Rec;
save ofdm_data BER_Rec;
save ofdm_data BER_Rec1;
%Watch the results
semilogy(EbN0_Rec,BER_Rec,'r-');
xlabel('Eb/N0 in dB');
ylabel('Bit error rate');
hold on
semilogy(EbN0_Rec,BER_Rec1,'*-');
xlabel('Eb/N0 in dB');
ylabel('hard decoding bit error rate');
hold on
legend('Zero forcing equalizing','MMSE equalizing');
grid