?? 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 + -