📄 vblast_ofdm.m
字号:
clear all;
close all;
fprintf('**** Simulation starts at: ');
disp(datestr(now,0));
% pdp : Multipath channel power delay profile.
% 相邻径之间的相对时延为一个采样周期
% L 多径数目: 0至L-1径。
L=3;
var0=1;
%%%% var0=(1-exp(-1))/(1-exp(-L)); % 多径总功率归一化
ll=0:L-1;
pdp=var0.*exp(-ll);
%%%%%%%%%%%%%%%%%%%%%%%%%%
% 在每个子载波的范围内分配能量与比特
%%%%%%%%%%%%%%%%%%%%%%%%%%%
num=4; % 平均每个天线的每个子载波分配的比特数
load parameters.txt; % 装入仿真参数
cases=parameters(1);
caseparas=parameters(2);
for kk=1:cases
% The number of transmit antennas.
Nt=parameters((kk-1)*caseparas+3);
% Nr: The number of receive antennas.
Nr=parameters((kk-1)*caseparas+4);
% N: The number of subcarriers.
N=parameters((kk-1)*caseparas+5);
% GI: CP length.
GI=parameters((kk-1)*caseparas+6);
% framenumber: The number of frames to simulate.
framenumber=parameters((kk-1)*caseparas+7);
% minSNR: Minimum SNR in dB to simulate begin with.
minSNR=parameters((kk-1)*caseparas+8);
% maxSNR: Maximum No in dB to simulate end with.
maxSNR=parameters((kk-1)*caseparas+9);
% stepSNR: SNR step in dB to increment.
stepSNR=parameters((kk-1)*caseparas+10);
% state: Modulation format: Adaptive=0; BPSK=1; QPSK=2; 16QAM=4; 64QAM=6.
state=parameters((kk-1)*caseparas+11);
% detection: Detection method: ZF=1; MMSE=2.
detection=parameters((kk-1)*caseparas+12);
% 线形星座预编码矩阵
% LCPindication=0; % 为零的时候不使用LCP预编码，用串行干扰删除；
% 若为5表示不用串行干扰删除，且不用LCP预编码。
% if LCPindication~=0
% switch LCPindication
% case 2
% lcp=[exp(-j*pi/4) exp(-j*5*pi/4)].';
% lcpint=[1 1].';
% LCP=1/sqrt(2)*[lcpint lcpint.*lcp];
% case 4
% lcp=[exp(-j*pi/8) exp(-j*5*pi/8) exp(-j*9*pi/8) exp(-j*13*pi/8)].';
% lcpint=[1 1 1 1].';
% LCP=1/2*[lcpint lcpint.*lcp lcpint.*(lcp.^2) lcpint.*(lcp.^3)];
% case 5
% LCP=eye(Nt);
% end
% invLCP=inv(LCP);
% end
fprintf('==== Case %d ====\n',kk);
disp('---- Simulation parameter list:');
fprintf('Transmit antennas: %d.\n', Nt);
fprintf('Receive antennas: %d.\n', Nr);
fprintf('The number of subcarriers : %d.\n', N);
fprintf('CP length : %d.\n', GI);
fprintf('Total frames: %d.\n', framenumber);
fprintf('SNR range to simulate: [%d:%1.1f:%d] dB.\n', minSNR,stepSNR,maxSNR);
switch(state)
case 0
disp('Adaptive modulation.');
fmod='ADPM';
case 1
disp('BPSK modulation.');
fmod='BPSK';
case 2
disp('QPSK modulation.');
fmod='QPSK';
case 4
disp('16QAM modulation.');
fmod='QAM16';
case 6
disp('64QAM modulation.');
fmod='QAM64';
end
switch detection
case 1
disp('ZF detection.');
case 2
disp('MMSE detection.');
end
% fprintf('LCPindication: %d.\n', LCPindication);
disp('----------------------------');
SNR=[minSNR:stepSNR:maxSNR];
randn('state',0);
rand('state',1);
frame_length=N+GI; % 每帧的长度;
err=zeros(1,length(SNR));
averBER=zeros(1,length(SNR));
Ptotmax=Nt*N;
if state==0
bitstotal=num*Nt*N*framenumber; % 每个发天线的每个子载波上平均分num个比特。对应比较2.^num阶调制。
else
bitstotal=Nt*N*framenumber*state; % 发送的总比特数
end
for SNRIndex=1:length(SNR)
sig2=Nt/(10.^(SNR(SNRIndex)/10)); % 此时的噪声功率为按整个ofdm符号平均能量归一化时的
% 每个发送天线每次发送一帧
last_symbol=zeros(Nt,length(pdp)-1); %信道记忆产生的影响
for kk=1:framenumber
% 完成一个帧的发送与接收
% 假设信道为准静态信道，在一个ofdm符号内保持不变
% 在符号间隔处发生跃变
[H H_f]=create_channel(Nt,Nr,pdp,N,GI);
% H_f为每个子载波对应的MIMO信道矩阵
% if state==0
% for ii=1:N
% Hnorm2(ii)=1/Nt/Nr*trace(H_f{ii}*H_f{ii}');
% end
% [Psubc,Rsubc]=BPallocRNew(Hnorm2/N/sig2,N,Ptotmax,N*num); % 每个发天线的每个子载波上平均分配的功率与比特数
% end
w=cell(1,N);
korder=cell(1,N);
% G2=cell(1,N); % 有LCP编码的时候用
if detection==2
eyenew=eye(Nr); % Detection method: ZF=1; MMSE=2.
end
for f=1:N
Hnew=H_f{f};
% if LCPindication
% if detection==1
% ZF
% G2{f}=pinv(Hnew); %pinv is the pseudoinverse
% else
% MMSE
% G2{f}=Hnew'*pinv(Hnew*Hnew'+N*sig2*eyenew);
% end
% else
korder{f}=zeros(1,Nt);
krest=[1:Nt];
for nt=1:Nt
if detection==1
% ZF
G=pinv(Hnew); %pinv is the pseudoinverse
if state==0
[k,postSNR]=argmaxzf(G,N*sig2,krest); % find the maximum norm form krest columes of G
else
[k,postSNR]=argminzf(G,N*sig2,krest); % find the minimum norm form krest columes of G
end
else
% MMSE
G=Hnew'*pinv(Hnew*Hnew'+N*sig2*eyenew);
if state==0
[k,postSNR]=argmaxmmse(G,Hnew,N*sig2,krest);
else
[k,postSNR]=argminmmse(G,Hnew,N*sig2,krest);
end
eyenew(:,k)=zeros(Nr,1);
end
% 查找出最小的SNR所在的列,用零替代,产生新的信道伪逆
korder{f}(nt)=k;
w{f}(k,:)=G(k,:);
krest=setdiff(krest,k); % delete the selected k
Hnew(:,k)=zeros(Nr,1);
if state==0
% averpostSNR(SNRIndex,f)=averpostSNR(SNRIndex,f)+postSNR;
SNRsubch(k,f)=postSNR;
end
end
end
% 产生发送数据
tr_data=cell(Nt,N);
symbol=zeros(Nt,N);
if state==0
SNRsubch2=reshape(SNRsubch.',1,Nt*N); % 变形后每个发天线的子载波信道对应信噪比占连续的一段。
[P1,R1]=BPallocRNew(SNRsubch2,Nt*N,Ptotmax,num*Nt*N);
P=reshape(P1,N,Nt).'; % 每个发天线的各子载波的能量占一行
R=reshape(R1,N,Nt).'; % 每个发天线的各子载波的比特数占一行
for nt=1:Nt
for f=1:N
tr_data{nt,f}=randn(1,R(nt,f))>0;
switch R(nt,f)
case 0
symbol(nt,f)=zeros(1,1);
case 1 % 根据自适应算法比特增量为2，此调制方式没有被使用
symbol(nt,f)=BPSKMod(tr_data{nt,f});
case 2
symbol(nt,f)=QPSKMod(tr_data{nt,f});
case 4
symbol(nt,f)=QAM16Mod(tr_data{nt,f});
case 6
symbol(nt,f)=QAM64Mod(tr_data{nt,f});
end
end
end
Proot=sqrt(P);
symbol=Proot.*symbol;
else
tr_data2=randn(Nt,N*state)>0;
for nt=1:Nt
symbol(nt,:)=modulate(tr_data2(nt,:),state);
end
end
% if LCPindication
% symbol=LCP*symbol;
% end
% ofdm调制，进行ifft变换
ofdm_symbol=ifft(symbol,[],2);
% 添加CP
tr_symbol=[ofdm_symbol(:,end-GI+1:end) ofdm_symbol];
% 发送信号平均功率归一化
tr_symbol = tr_symbol*N/sqrt(N);
% 通过信道
re_cp_symbol= channel(H,tr_symbol,sig2,last_symbol,Nr,Nt,pdp);
last_symbol=tr_symbol(:,end-length(pdp)+2:end);
% 恢复信号的功率水平
re_cp_symbol=re_cp_symbol/sqrt(N);
% 去除cp
re_symbol=re_cp_symbol(:,GI+1:end);
% ofdm解调，进行fft
de_ofdm=fft(re_symbol,[],2);
% 对每个载波进行vblast检测
bitRx=cell(Nt,N); % 存放自适应的接收比特
y=zeros(Nt,N);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% if LCPindication
% for f=1:N
% y(:,f)=invLCP*G2{f}*de_ofdm(:,f);
% for nt=1:Nt
% if state==0
% switch Rsubc(f)
% case 0
% case 2
% bitRx{nt,f}=QPSKDemod(y(nt,f));
% case 4
% bitRx{nt,f}=QAM16Demod(y(nt,f));
% case 6
% bitRx{nt,f}=QAM64Demod(y(nt,f));
% end
%else
% bitRx{nt,f}=demodulate(y(nt,f),state);
%end
%end
%end
%end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for f=1:N
symbolRx=de_ofdm(:,f);
for nt=1:Nt
y=w{f}(korder{f}(nt),:)*symbolRx;
% decision and demodulation
if state==0
switch R(korder{f}(nt),f)
case 0
a=0;
case 2
bitRx{korder{f}(nt),f}=QPSKDemod(y/Proot(korder{f}(nt),f));
a=QPSKMod(bitRx{korder{f}(nt),f})*Proot(korder{f}(nt),f);
case 4
bitRx{korder{f}(nt),f}=QAM16Demod(y/Proot(korder{f}(nt),f));
a=QAM16Mod(bitRx{korder{f}(nt),f})*Proot(korder{f}(nt),f);
case 6
bitRx{korder{f}(nt),f}=QAM64Demod(y/Proot(korder{f}(nt),f));
a=QAM64Mod(bitRx{korder{f}(nt),f})*Proot(korder{f}(nt),f);
end
else
bitRx{korder{f}(nt),f}=demodulate(y,state);
a=modulate(bitRx{korder{f}(nt),f},state);
end
symbolRx=symbolRx-H_f{f}(:,korder{f}(nt))*a;
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
tr=[];
if state==0
for kk=1:Nt
for f=1:N
tr=[tr tr_data{kk,f}];
end
end
else
tr=reshape(tr_data2.',1,Nt*N*state);
end
Rx=[];
for kk=1:Nt
for f=1:N
Rx=[Rx bitRx{kk,f}];
end
end
err(SNRIndex)=err(SNRIndex)+sum(tr~=Rx); % 加上一帧中误比特数
end % 对应一帧检测完
end % 对应SNRIndex所有帧都检测完
averBER=err/bitstotal;
if state==0
adap4BER11=averBER;
save 'adap4BER11.mat' SNR adap4BER11;
else
qam16BER11=averBER;
save 'qam16BER11.mat' SNR qam16BER11;
end
end % End one case
% switch LCPindication
% case 2
% lcp22BER=averBER;
% save 'lcp22BER.mat' SNR averBER
% case 4
% lcp44BER=averBER;
% save 'lcp44BER.mat' SNR averBER
% case 5
% nolcp_canBER44=averBER;
% save 'nolcp_canBER44' SNR averBER
% end
figure;
semilogy(SNR,qam16BER11,'bo--',SNR,adap4BER11,'r^-');
axis([0 30 10.^-4 1]);
xlabel('average SNR (dB)');
ylabel('Raw BER');
legend('conventional','adptive');
grid on;
fprintf('\n**** Ending time: ');
disp(datestr(now,0));
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -