MIMO_OFDM.m
上传用户:m_sun_001
上传日期:2014-07-30
资源大小:1115k
文件大小:6k
源码类别:

matlab例程

开发平台:

Matlab

  1. %------------------------------------------
  2. % EE359 final project, Fall 2002
  3. % Channel estimation for a MIMO-OFDM system
  4. % By Shahriyar Matloub               
  5. %------------------------------------------
  6. clear all;
  7. %close all;
  8. i=sqrt(-1);
  9. Rayleigh=1;
  10. AWGN=0;                             % for AWGN channel 
  11. MMSE=0;                             % estimation technique
  12. Nsc=64;                             % Number of subcarriers
  13. Ng=16;                              % Cyclic prefix length
  14. SNR_dB=[0 5 10 15 20 25 30 35 40];  % Signal to noise ratio
  15. Mt=2;                               % Number of Tx antennas
  16. Mr=2;                               % Number of Rx antennas
  17. pilots=[1:Nsc/Ng:Nsc];              % pilot subcarriers 
  18. DS=5;                              % Delay spread of channel
  19. iteration_max=200;
  20. %%%%%%%%%%%%%%%%%%%%%%%%%%%%
  21. % Channel impulse response %
  22. %%%%%%%%%%%%%%%%%%%%%%%%%%%%
  23. if (Rayleigh)
  24.     N=50;
  25.     fm=100;
  26.     B=20e3;
  27.     fd=(rand(1,N)-0.5)*2*fm;
  28.     theta=randn(1,N)*2*pi;
  29.     c=randn(1,N);
  30.     c=c/sum(c.^2);
  31.     t=0:fm/B:10000*fm/B;
  32.     Tc=zeros(size(t));
  33.     
  34.     Ts=zeros(size(t));
  35.     for k=1:N
  36.        Tc=c(k)*cos(2*pi*fd(k)*t+theta(k))+Tc;
  37.        Ts=c(k)*sin(2*pi*fd(k)*t+theta(k))+Ts;
  38.     end
  39.     r=ones(Mt*Mr,1)*(Tc.^2+Ts.^2).^0.5;
  40.     index=floor(rand(Mt*Mr,DS)*5000+1);
  41. end
  42. MEE1=zeros(1,length(SNR_dB));
  43. MEE2=zeros(1,length(SNR_dB));
  44. for snrl=1:length(SNR_dB)
  45.     snrl
  46.     estimation_error1=zeros(Mt*Mr,Nsc);
  47.     estimation_error2=zeros(Mt*Mr,Nsc);
  48.     R1=besselj(0,2*pi*fm*(Nsc+Ng)/B);
  49.     sigma2=10^(-SNR_dB(snrl)/10);
  50.     aa=(1-R1^2)/(1-R1^2+sigma2);
  51.     bb=sigma2*R1/(1-R1^2+sigma2);
  52.     for iteration=1:iteration_max
  53.         %iteration    
  54.         if AWGN==1
  55.             h=ones(Mt*Mr,1);
  56.         else
  57.             phi=rand*2*pi;
  58.             h=r(index+iteration)*exp(j*phi);
  59.             %h=rand(Mt*Mr,DS);
  60.             h=h.*(ones(Mt*Mr,1)*(exp(-0.5).^[1:DS]));
  61.             h=h./(sqrt(sum(abs(h).^2,2))*ones(1,DS));
  62.         end
  63.         CL=size(h,2);                                               % channel length
  64.         data_time=zeros(Mt,Nsc+Ng);
  65.         data_qam=zeros(Mt,Nsc);
  66.         data_out=zeros(Mr,Nsc);
  67.         output=zeros(Mr,Nsc);
  68.         for tx=1:Mt
  69.             data_b=0*round(rand(4,Nsc));                                  % data
  70.             data_qam(tx,:)=j*(2*(mod(data_b(1,:)+data_b(2,:),2)+2*data_b(1,:))-3)+...
  71.             2*(mod(data_b(3,:)+data_b(4,:),2)+2*data_b(3,:))-3;
  72.             for loop=1:Mt 
  73.                 data_qam(tx,pilots+loop-1)=(1+j)*(loop==tx);              % pilots
  74.             end
  75.             data_time_temp=ifft(data_qam(tx,:));
  76.             data_time(tx,:)=[data_time_temp(end-Ng+1:end) data_time_temp];
  77.         end
  78.     
  79.         for rx=1:Mr
  80.             for tx=1:Mt
  81.                 output_temp=conv(data_time(tx,:),h((rx-1)*Mt+tx,:));
  82.                 output(rx,:)=output_temp(Ng+1:Ng+Nsc)+output(rx,:);
  83.             end
  84.             np=(sum(abs(output(rx,:)).^2)/length(output(rx,:)))*sigma2;
  85.             noise=(randn(size(output(rx,:)))+i*randn(size(output(rx,:))))*sqrt(np);
  86.             output(rx,:)=output(rx,:)+noise;
  87.             data_out(rx,:)=fft(output(rx,:));
  88.         end
  89. %%%%%%%%%%%%%%%%%%%%%%
  90. % Channel estimation %
  91. %%%%%%%%%%%%%%%%%%%%%%
  92.     
  93.         H_act=zeros(Mt*Mr,Nsc);
  94.         H_est1=zeros(Mt*Mr,Nsc);
  95.         H_est2=zeros(Mt*Mr,Nsc);
  96.         i=1;
  97.         for tx=1:Mt
  98.             for rx=1:Mr
  99.                 H_est_temp=data_out(rx,pilots+tx-1)./data_qam(tx,pilots+tx-1);
  100.                 %H_est_temp2=aa*abs(H_est_temp1)+bb*abs(H_est2((rx-1)*Mt+tx,:));
  101.                 h_time=ifft(H_est_temp);
  102.                 h_time=[h_time zeros(1,Nsc-length(h_time))];
  103.                 H_est1((rx-1)*Mt+tx,:)=fft(h_time);
  104.                 H_est2((rx-1)*Mt+tx,:)=((aa*abs(H_est1((rx-1)*Mt+tx,:))+bb*abs(H_est2((rx-1)*Mt+tx,:)))...
  105.                     .*H_est1((rx-1)*Mt+tx,:))./abs(H_est1((rx-1)*Mt+tx,:));
  106.                 if (tx>1)
  107.                     H_est1((rx-1)*Mt+tx,:)=[H_est1((rx-1)*Mt+tx,Nsc-tx+2:Nsc) H_est1((rx-1)*Mt+tx,1:Nsc-tx+1)];
  108.                     H_est2((rx-1)*Mt+tx,:)=[H_est2((rx-1)*Mt+tx,Nsc-tx+2:Nsc) H_est2((rx-1)*Mt+tx,1:Nsc-tx+1)];    
  109.                 end
  110.                 H_act((rx-1)*Mt+tx,:)=fft([h((rx-1)*Mt+tx,:) zeros(1,Nsc-CL)]);
  111.                 error1=(abs(H_act((rx-1)*Mt+tx,:)-H_est1((rx-1)*Mt+tx,:)).^2);
  112.                 error2=(abs(H_act((rx-1)*Mt+tx,:)-H_est2((rx-1)*Mt+tx,:)).^2);
  113.                 %error=(abs(H_act((rx-1)*Mt+tx,:)-H_est((rx-1)*Mt+tx,:)).^2)./(abs(H_act((rx-1)*Mt+tx,:)).^2);
  114.                 estimation_error1((rx-1)*Mt+tx,:)=estimation_error1((rx-1)*Mt+tx,:)+error1;                 
  115.                 estimation_error2((rx-1)*Mt+tx,:)=estimation_error2((rx-1)*Mt+tx,:)+error2; 
  116.                 %subplot(Mt*Mr,3,i),plot([0:Nsc-1],abs(H_act((rx-1)*Mt+tx,:))); i=i+1;
  117.                 %subplot(Mt*Mr,3,i),plot([0:Nsc-1],abs(H_est((rx-1)*Mt+tx,:))); i=i+1;
  118.                 %subplot(Mt*Mr,3,i),plot([0:Nsc-1],abs(error)); i=i+1;
  119.             end
  120.         end  
  121.     end
  122.     estimation_error1=estimation_error1/iteration_max;
  123.     estimation_error2=estimation_error2/iteration_max;
  124.     %estimation_error=min(estimation_error,10*iteration_max*ones(size(estimation_error)));
  125.     %for i=1:Mt*Mr
  126.     %    subplot(Mt*Mr,2,2*i-1),plot([0:Nsc-1],estimation_error1(i,:));    
  127.     %    subplot(Mt*Mr,2,2*i),plot([0:Nsc-1],estimation_error2(i,:));
  128.     %end
  129.     MEE1(snrl)=sum(sum(estimation_error1))/(Mt*Mr*Nsc);
  130.     MEE2(snrl)=sum(sum(estimation_error2))/(Mt*Mr*Nsc);
  131. end
  132. plot(SNR_dB,10*log10(MEE1));    
  133. hold on;
  134. plot(SNR_dB,10*log10(MEE2),'r');
  135. %H_act=fft([h_zeros(1,Nsc-CL)]).';
  136. error1=(abs(H_act-H_est1).^2)./(abs(H_act).^2);
  137. error2=(abs(H_act-H_est2).^2)./(abs(H_act).^2);
  138. %%%%%%%%%
  139. % Plots %
  140. %%%%%%%%%
  141. fig=4;
  142. i=1;
  143. subplot(fig,1,i),plot([0:length(H_act)-1],abs(H_act));    i=i+1;
  144. subplot(fig,1,i),plot([0:length(H_est1)-1],abs(H_est1));  i=i+1;
  145. subplot(fig,1,i),plot([0:length(H_est2)-1],abs(H_est2));  i=i+1;
  146. subplot(fig,1,i),plot([0:length(error1)-1],error1);       i=i+1;
  147. subplot(fig,1,i),plot([0:length(error2)-1],error2);