MATLAB SYSTEM VIEW通信原理实验及系统仿真

clear all close all clc i=sqrt(-1); CarryN=100; %子载波数 IFFTLen=512; %傅立叶变换长度 SymbN=66; %符号数/载波 CPLen=10; %循环前缀长度 AddPreLen=IFFTLen+CPLen;%添加循环前缀后的长度 Mpsk=4; %QPSK基带调制方式， BitPerSym=log2(Mpsk); %比特数/符号 % 调制符号数为4，发射天线数为3的4*3复正交空时编码矩阵[x1 -x2 -x3;x2* x1* 0;x3* 0 x1*;0 -x3* x2*] SendMetr=[1 -2 -3;2+j 1+j 0;3+j 0 1+j;0 -3+j 2+j]; co_time=size(SendMetr,1); %发射一组符号的时间周期：4 Nt=size(SendMetr,2); %发射天线数目：3 Nr=2; %接收天线数目：2 snr_min=3; %最小信噪比 ，dB snr_max=15; %最大信噪比 ，dB graph_inf_bit=zeros(snr_max-snr_min+1,2,Nr); %绘图信息存储矩阵 graph_inf_sym=zeros(snr_max-snr_min+1,2,Nr); num_X=1; for cc_ro=1:co_time%找出编码矩阵中最大的调制符号的序号 for cc_co=1:Nt num_X=max(num_X,abs(real(SendMetr(cc_ro,cc_co)))); end end co_x=zeros(num_X,1);%用于统计各个调制符号的个数 for con_ro=1:co_time for con_co=1:Nt%用于确定矩阵“SendMetr”中元素的位置，符号以及共轭情况 if abs(real(SendMetr(con_ro,con_co)))~=0 delta(con_ro,abs(real(SendMetr(con_ro,con_co))))=sign(real(SendMetr(con_ro,con_co))); %判断各个编码符号的正负 epsilon(con_ro,abs(real(SendMetr(con_ro,con_co))))=con_co;%天线号 co_x(abs(real(SendMetr(con_ro,con_co))),1)=co_x(abs(real(SendMetr(con_ro,con_co))),1)+1; %相同调制符号的个数 eta(abs(real(SendMetr(con_ro,con_co))),co_x(abs(real(SendMetr(con_ro,con_co))),1))=con_ro;%各符号传输时隙 coj_mt(con_ro,abs(real(SendMetr(con_ro,con_co))))=imag(SendMetr(con_ro,con_co)); %虚部 end end end eta=eta.'; eta=sort(eta); %Sort is not necessary eta=eta.'; carriers = (1: CarryN) + (floor(IFFTLen/4) - floor(CarryN/2));%载波坐标 79：178 conjugate_carriers=IFFTLen-carriers+2;%载波变换坐标 435：336 tx_training_symbols=training_symbol(Nt,CarryN);%准备训练序列 baseband_out_length = CarryN * SymbN; %基带输出长度 %Start simulation snrlength=snr_max-snr_min; x_sym=zeros(snrlength,1); y_sym=zeros(snrlength,1); n_sym=zeros(snrlength,1); m_sym=zeros(snrlength,1); for itNum=1: 10 %迭代仿真次数：10 for SNR=snr_min:snr_max n_err_sym=zeros(1,Nr); n_err_bit=zeros(1,Nr); Perr_sym=zeros(1,Nr); Perr_bit=zeros(1,Nr); re_met_sym_buf=zeros(CarryN,SymbN,Nr); re_met_bit=zeros(baseband_out_length,BitPerSym,Nr); baseband_out=round(rand(baseband_out_length,BitPerSym)); %生成随机输入序列 de_data=bi2de(baseband_out); %二进制向十进制转换 data_buf=pskmod(de_data,Mpsk,0); %QPSK调制 ，Mpsk=4 carrier_matrix=reshape(data_buf,CarryN,SymbN); %按子载波数进行串并转换 %取数为空时编码做准备，此处每次取每个子载波上连续的两个数 for tt=1:Nt:SymbN data=[]; for ii=1:Nt tx_buf_buf=carrier_matrix(:,tt+ii-1); data=[data;tx_buf_buf]; end XX=zeros(co_time*CarryN,Nt);%用于存放空时编码后的符号 for con_r=1:co_time %本循环进行空时编码 for con_c=1:Nt if abs(real(SendMetr(con_r,con_c)))~=0 if imag(SendMetr(con_r,con_c))==0 XX((con_r-1)*CarryN+1:con_r*CarryN,con_c)=data((abs(real(SendMetr(con_r,con_c)))-1)*CarryN+1:abs(real(SendMetr(con_r,con_c)))... *CarryN,1)*sign(real(SendMetr(con_r,con_c))); else %如果虚部不为0，求复共轭 XX((con_r-1)*CarryN+1:con_r*CarryN,con_c)=conj(data((abs(real(SendMetr(con_r,con_c)))-1)*CarryN+1:abs(real(SendMetr(con_r,con_c)))... *CarryN,1))*sign(real(SendMetr(con_r,con_c))); end end end end %空时编码结束 XX=[tx_training_symbols;XX]; %插入导频，即训练序列 rx_buf=zeros(1,AddPreLen*(co_time+1),Nr); %准备各天线接收缓冲 for rev=1:Nr for ii=1:Nt tx_buf=reshape(XX(:,ii),CarryN,co_time+1); %各天线串并转换，准备IFFT IFFT_tx_buf=zeros(IFFTLen,co_time+1); IFFT_tx_buf(carriers,:)=tx_buf(1:CarryN,:); %准备IFFT复数序列 IFFT_tx_buf(conjugate_carriers,:)=conj(tx_buf(1:CarryN,:));%准备IFFT的共轭对称序列 time_matrix=ifft(IFFT_tx_buf); %IFFT运算 time_matrix=[time_matrix((IFFTLen-CPLen+1):IFFTLen,:);time_matrix];%插入循环前缀 tx=time_matrix(:)';%ii天线数据并串转换，准备发射 tx_tmp=tx; d=[4,5,6,2;4,5,6,2;4,5,6,2;4,5,6,2]; %模拟4径传输中各径时延 a=[0.2,0.3,0.4,0.5;0.2,0.3,0.4,0.5;0.2,0.3,0.4,0.5;0.2,0.3,0.4,0.5];%模拟4径传输中各径衰落 for jj=1:size(d,2) copy=zeros(size(tx)) ; for kk = 1 + d(ii,jj): length(tx) copy(kk) = a(ii,jj)*tx(kk - d(ii,jj)) ; end tx_tmp=tx_tmp+copy; %多径传输后的信号 end txch=awgn(tx_tmp,SNR,'measured'); %AWGN信道，添加高斯白噪声 rx_buf(1,:,rev)=rx_buf(1,:,rev)+txch; %接收数据 end rx_spectrum=reshape(rx_buf(1,:,rev),AddPreLen,co_time+1); %串并转换，准备各路FFT rx_spectrum=rx_spectrum(CPLen+1:AddPreLen,:); %去除循环前缀 FFT_tx_buf=zeros(IFFTLen,co_time+1); FFT_tx_buf=fft(rx_spectrum); %FFT运算 spectrum_matrix=FFT_tx_buf(carriers,:); %分离各子载波数据，含导频 Y_buf=(spectrum_matrix(:,2:co_time+1)); %去除导频 Y_buf=conj(Y_buf'); %求复共轭 spectrum_matrix1=spectrum_matrix(:,1); %取出导频 Wk=exp((-2*pi/CarryN)*i);%DFT基波 L=10; p=zeros(L*Nt,1); %进行基于导频的信道估计 for jj=1:Nt for l=0:L-1 for kk=0:CarryN-1 p(l+(jj-1)*L+1,1)=p(l+(jj-1)*L+1,1)+spectrum_matrix1(kk+1,1)*conj(tx_training_symbols(kk+1,jj))*Wk^(-(kk*l)); end end end h=p/CarryN; H_buf=zeros(CarryN,Nt); for ii=1:Nt for kk=0:CarryN-1 for l=0:L-1 H_buf(kk+1,ii)=H_buf(kk+1,ii)+h(l+(ii-1)*L+1,1)*Wk^(kk*l); end end end H_buf=conj(H_buf'); RRR=[]; for kk=1:CarryN %本循环完成空时解码 Y=Y_buf(:,kk); %接收到的数据 H=H_buf(:,kk); %估计的信道矩阵 for co_ii=1:num_X for co_tt=1:size(eta,2) if eta(co_ii,co_tt)~=0 if coj_mt(eta(co_ii,co_tt),co_ii)==0 %虚部为零 r_til(eta(co_ii,co_tt),:,co_ii)=Y(eta(co_ii,co_tt),:); a_til(eta(co_ii,co_tt),:,co_ii)=conj(H(epsilon(eta(co_ii,co_tt),co_ii),:)); else%虚部不为零，求复共轭 r_til(eta(co_ii,co_tt),:,co_ii)=conj(Y(eta(co_ii,co_tt),:)); a_til(eta(co_ii,co_tt),:,co_ii)=H(epsilon(eta(co_ii,co_tt),co_ii),:); end end end end RR=zeros(num_X,1); for iii=1:num_X %%本循环完成调制符号的估计 for ttt=1:size(eta,2) if eta(iii,ttt)~=0 RR(iii,1)=RR(iii,1)+r_til(eta(iii,ttt),1,iii)*a_til(eta(iii,ttt),1,iii)*delta(eta(iii,ttt