%******Pedestrian-A 信道BER仿真********%
clc
clear all
%% 参数设置
N=128; %子载波数目
M=5; %每个子载波里的字符号数目
qam_M=64; %QAM调制阶数
psk_M=4; %PSK调制阶数
mod_M=[qam_M psk_M];
g_ori=zeros(2,M*N+1);
g_ori_final=zeros(2,M*N);
B_all=zeros(2*M*N,2*M*N);
B=mat2cell(B_all,[M*N,M*N],2*M*N);
%% 建立PA信道
ts=1/640;
fd=8.3*2/3.8;
tau=[0 200 800 1200 2300 3700]*10^(-9);
pdb=[0 -0.9 -4.9 -8.0 -7.8 -23.9];
chan=comm.RayleighChannel('SampleRate',1/ts, ...
'PathDelays',tau, ...
'AveragePathGains',pdb, ...
'MaximumDopplerShift',fd);
%%
%% 建立滤波器
beta=[0.5 0.9]; % 滚降系数
for ig=1:2
g_ori(ig,:)=rcosdesign(beta(ig),M,N); % g_ori为升余弦滚降滤波器（1×641）
g_ori_final(ig,:)=g_ori(ig,2:641);
g(:,ig)=g_ori_final(ig,:)'; % 将g_ori转置
G=zeros(M*N,2*M); % 构造矩阵G
for i=1:10 % 对g进行循环移位
G(:,i)=circshift(g(:,ig),64*(i-1));
end
[V,L,U]=svd(G); % L对角元素是对G'G进行特征分解得到特征值矩阵对角元素的平方根，特征值降序排列
L_diag=diag(L);
L_unnegative=L_diag(1:5,:); % 取L的前五个非0值
P = zeros(10,10); %建立功率分配矩阵P
P(1:5,1:5)=diag(1./L_unnegative);
P(6:10,6:10)=0; % P(2M×2M)(10×10)
Qpre=U; %因为G特征分解得到的矩阵U应与G'G特征分解得到的特征向量矩阵相等，故直接取Qpre=U
Qpost=U';
I=eye(128,128);
P_bl=kron(I,P);
U_bl=kron(I,U);
A=zeros(M*N,2*M*N);
w=zeros(N,1);
for iw=1:128
for jw=0:127
w(jw+1)=exp(1j*2*pi*jw*iw/N);
end
A(:,10*(iw-1)+1:10*iw)=diag([w' w' w' w' w']')*G;
end
B(ig)={A*U_bl*P_bl};
end
B1=B(1); %滚降系数为0.5的GFDM调制矩阵
B2=B(2); %滚降系数为0.9的GFDM调制矩阵
%%
times=10;
ebno=0:10;
num_err_all=zeros(2*times,2*length(ebno));
num_err=mat2cell(num_err_all,[times times],[length(ebno) length(ebno)]);
E2R=[10*log10((6)) 10*log10((2))]; %Eb/N0与SNR之间的转换因子
for ix=1:times
%*****产生信号源******%
x_ori_qam=randi([0, 1], 5*log2(qam_M), 128); % 随机生成一个0、1组成的30行128列数据矩阵x_ori(M*Mlog2(QAM)×N)(30×128)
x_con_qam=qammod(x_ori_qam, qam_M, 'bin', 'InputType', 'bit', 'OutputDataType', numerictype(1,16,10)); % 对x_ori进行64QAM调制得到信号x_con(M×N)(5×128)
x_ori_psk=randi([0, 1], 5*log2(psk_M), 128); % 随机生成一个0、1组成的10行128列数据矩阵x_ori(M*Mlog2(QAM)×N)(10×128)
x_con_psk=qammod(x_ori_psk, psk_M, 'bin', 'InputType', 'bit', 'OutputDataType', numerictype(1,16,10)); % 对x_ori进行64QAM调制得到信号x_con(M×N)(5×128)
x_ori_all=[x_ori_qam;x_ori_psk];
x_ori=mat2cell(x_ori_all,[30, 10]);
%*****对x进行补零*****
o=zeros(M,N);
D_qam=[x_con_qam;o]; % x(2M×N)(10×128)
d_qam=reshape(D_qam,2*M*N,1); % d(1280×1)
D_psk=[x_con_psk;o]; % x(2M×N)(10×128)
d_psk=reshape(D_psk,2*M*N,1); % d(1280×1)
s_all=zeros(1280,2);
s=mat2cell(s_all,[640 640],[1 1]);
s(1,1)={cell2mat(B1)*d_qam}; %QAM with 滚降系数0.5
s(1,2)={cell2mat(B2)*d_qam}; %QAM with 滚降系数0.9
s(2,1)={cell2mat(B1)*d_psk}; %PSK with 滚降系数0.5
s(2,2)={cell2mat(B2)*d_psk}; %PSK with 滚降系数0.9
% 通过PA信道 并加噪声
F=zeros(1,length(ebno));
K=zeros(1,length(ebno));
for isr=1:2
for isc=1:2
ss=cell2mat(s(isr,isc));
for i=1:length(ebno)
snr= E2R(isr)+(i-1);
sigPower=10*log10(sum(abs(ss.data).^2)/length(ss.data));
s_pa=chan(double(ss)); %通过PA信道
H=fft(s_pa)./fft(ss.data);
s_add_noise=awgn(s_pa,snr,sigPower); %加高斯白噪声
eq_signal=fft(s_add_noise)./H;
s_recovered=cell2mat(B(isc))'*ifft(eq_signal);
d_recovered=reshape(s_recovered,2*M,N);
x_recovered=d_recovered(1:M,:); %去零
x_demod=qamdemod(x_recovered,mod_M(isr), 'bin', 'OutputType', 'bit');
[F(i),K(i)]=symerr(cell2mat(x_ori(isr)),x_demod);
end
num_err{isr,isc}(ix,:)=F;
end
end
end
figure(1)
ber_all=zeros(2,2*length(ebno));
ber=mat2cell(ber_all,[1 1],[length(ebno) length(ebno)]);
ber(1,1)={sum(cell2mat(num_err(1,1)))/(times*5*log2(qam_M)*N)};
ber(1,2)={sum(cell2mat(num_err(1,2)))/(times*5*log2(qam_M)*N)};
ber(2,1)={sum(cell2mat(num_err(2,1)))/(times*5*log2(psk_M)*N)};
ber(2,2)={sum(cell2mat(num_err(2,2)))/(times*5*log2(psk_M)*N)};
semilogy(ebno,cell2mat(ber(1,1)),'-r');
hold on
semilogy(ebno,cell2mat(ber(1,2)),'or');
hold on
semilogy(ebno,cell2mat(ber(2,1)),'--b');
hold on
semilogy(ebno,cell2mat(ber(2,2)),'ob');
legend('β=0.5,64QAM','β=0.9,64QAM','β=0.5,QPSK','β=0.9,QPSK');
xlabel('Eb/No [dB]'); ylabel('BER');
title('Pedestrian-A channel');
grid();
axis([0 10 10^(-6) 1])
