tic
clear all;
clc;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
%Initial Configuration%
QAMOrder=4;
N=64;
Lcp=8;
iter_num = 2;
% CPsize=0;
BandWidt=30720000;
Simu_Num = 10e6;
sample_time=1/BandWidt;%sample_time=1e-7;
DataPerFrame=400;
loop_length=ceil(Simu_Num/(N*DataPerFrame));%ceil((BitCalcNumb/log2(QAMOrder))/FFTsize);
h = modem.qammod('M', QAMOrder, 'PhaseOffset', 0, 'SymbolOrder',...
'Gray', 'InputType', 'Bit');
hDemo=modem.qamdemod('M',QAMOrder,'SymbolOrder','Gray','OutputType','Bit');
%信道参数设置。
% %Define the channel
Path_len=32;
Path_Gain=zeros(1,Path_len);
for n=1:16
Path_Gain(n)=sqrt(0.70/16);
end
for n=16+1:Path_len
Path_Gain(n)=sqrt(0.30/(Path_len-16));
end
%main program
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
EbN0 = 4:4:16;
BER=zeros(1,length(EbN0));
BER_CPR=zeros(1,length(EbN0),iter_num+1);
for EbN0_index=1:length(EbN0)
bit_err_sum=0;
bit_CPR_err_sum=zeros(1,iter_num+1);
for loop=1:loop_length
%初始bit数据data
Tx_bit_for_comp=randi([0,1],[N*log2(QAMOrder),DataPerFrame]);
%QAM调制。
Data_Tx_FD=modulate(h,Tx_bit_for_comp);
% scatterplot(reshape(temp1,1,numel(temp1)));
%ifft.
Data_Tx_TD =(sqrt(N)).*ifft(Data_Tx_FD,N,1);
%加CP.
Data_Tx_TD_CP = zeros(N+Lcp,DataPerFrame);
Data_Tx_TD_CP(1:Lcp,:) = Data_Tx_TD(end-Lcp+1:end,:);
Data_Tx_TD_CP(Lcp+1:end,:) = Data_Tx_TD;
%并串转换。
OutDataPoweEsti=reshape(Data_Tx_TD_CP,1,numel(Data_Tx_TD_CP));
%%
%数据过信道
%%pass fading channel
L1=Path_len;
ch=(randn(1,L1)+1i*randn(1,L1)).*Path_Gain/sqrt(2); %原始信道系数
signal_fading_serial_col = conv(ch.',OutDataPoweEsti);
Data_trans_multi_chan = signal_fading_serial_col(1:length(OutDataPoweEsti)); %过信道后进行截断
%add AWGN
y = OutDataPoweEsti; %这样处理是否可以???
[M1,N1]=size(y);
x_energ=sum(sum(abs(y).^2))/M1/N1;
Eb = x_energ/log2(QAMOrder);
n_energ=10^(-EbN0(EbN0_index)/10)*Eb; %对数EbN0到线性坐标的变换
n2=(randn(M1,N1)+1i*randn(M1,N1))*sqrt(n_energ/2); %噪声
Tmp_Nv=sum(sum(abs(n2.*n2)))/N1; % Estimate of noise variance
signal_awgn = Data_trans_multi_chan+n2;
% signal_awgn = Data_trans_multi_chan;%不加噪声;
%% 接收端处理
Data_Rx_serial=signal_awgn;
Data_Rx = reshape(Data_Rx_serial,N+Lcp,DataPerFrame);%串并转换
%去CP.
Data_Rx_CPremoved = Data_Rx(Lcp+1:end,:);
% 常规解调 ******************************************************************************************
%fft
Data_Rx_FD = (1/sqrt(N)).*fft(Data_Rx_CPremoved, N, 1);
%均衡.
Data_Rx_FD_EQ = zeros(N,DataPerFrame);
CH = fft(ch.', N, 1);
for count_col=1:DataPerFrame
Data_Rx_FD_EQ(:,count_col)=Data_Rx_FD(:,count_col)./CH ;
end
% tmpinData=reshape(Data_Rx_FD_EQ,N*DataPerFrame,1);
% Datascatter=scatterplot(tmpinData);
Rx_bit_for_comp=demodulate(hDemo,Data_Rx_FD_EQ);
[number1,ratio1]=biterr(Tx_bit_for_comp,Rx_bit_for_comp);
bit_err_sum=bit_err_sum+ratio1;
% CPR处理 ******************************************************************************************
L=Path_len-1;
% Xj = Data_Rx_FD_EQ;
% xj = (sqrt(N)).*ifft(Xj,N,1);
xj = zeros(N,DataPerFrame,iter_num+1);
for i_iter = 1:iter_num+1
xj(:,1,i_iter) = Data_Tx_TD(:,1); %第一个符号假设已知
end
r = Data_Rx_CPremoved;
r1 = Data_Rx; %注意到r1是未去除CP的接收信号
G=Lcp;
W_factor = zeros(L-G,1);
h_pow = ch.*(ch').';
h_pow_sum = sum(h_pow);
for k=0:L-G-1
W_factor(k+1) = sum(h_pow(G+k+2:L+1))/h_pow_sum; %文中公式(7)
end
for ii_iter = 1:iter_num+1
if ii_iter == 1
xj(:,:,1) = CPR_initialize(h,hDemo,N,DataPerFrame,ch,xj,r,r1,W_factor,L,G); %计算出第一次解调的数据
else
xj(:,:,ii_iter) = CPR_iteration(h,hDemo,N,DataPerFrame,ch,xj,r,L,G,ii_iter); %计算出第n次迭代后解调的数据
end
%fft
Data_CPR_Rx_TD = xj(:,:,ii_iter);
Data_CPR_Rx_FD = (1/sqrt(N)).*fft(Data_CPR_Rx_TD, N, 1);
Rx_CPR_bit_for_comp=demodulate(hDemo,Data_CPR_Rx_FD);
[number_CPR,ratio_CPR]=biterr(Tx_bit_for_comp(:,2:end-1),Rx_CPR_bit_for_comp(:,2:end-1));
bit_CPR_err_sum(ii_iter)=bit_CPR_err_sum(ii_iter)+ratio_CPR;
end
end
BER(1,EbN0_index)=bit_err_sum/loop_length;
for iii_iter = 1:iter_num+1
BER_CPR(1,EbN0_index,iii_iter)=bit_CPR_err_sum(iii_iter)/loop_length;
end
display(['EbN0_in=',num2str(EbN0(EbN0_index))]);
end
figure(1);
semilogy(EbN0,BER,'-k+',EbN0,BER_CPR(:,:,1),'-ro',EbN0,BER_CPR(:,:,2),'-b*',EbN0,BER_CPR(:,:,3),'-g.');
box('on');grid('on');hold('all');
legend('No process','CPR I=0','CPR I=1','CPR I=2');
xlabel('EbN0(dB)');ylabel('BER');
title(['PathLen=',num2str(Path_len),',','CPLen=',num2str(Lcp),',','fftsize=',num2str(N),',','frame=',num2str(DataPerFrame),',','QAMorder=',num2str(QAMOrder),',','SimuNum=',num2str(Simu_Num)]);
toc
%
%%%%%%%%%%%%Simulting Loop%%%%%%%%%%%%%%%%%%%%%%%%%
%Programmming end.