%------------------------------------------
%?EE359?final?project?Fall?2002
%?Channel?estimation?for?a?MIMO-OFDM?system
%?By?Shahriyar?Matloub???????????????
%------------------------------------------

clear?all;
%close?all;
i=sqrt（-1）;
Rayleigh=1;
AWGN=0;?????????????????????????????%?for?AWGN?channel?
MMSE=0;?????????????????????????????%?estimation?technique
Nsc=64;?????????????????????????????%?Number?of?subcarriers
Ng=16;??????????????????????????????%?Cyclic?prefix?length
SNR_dB=[0?5?10?15?20?25?30?35?40];??%?Signal?to?noise?ratio
Mt=2;???????????????????????????????%?Number?of?Tx?antennas
Mr=2;???????????????????????????????%?Number?of?Rx?antennas
pilots=[1:Nsc/Ng:Nsc];??????????????%?pilot?subcarriers?
DS=5;??????????????????????????????%?Delay?spread?of?channel
iteration_max=200;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%?Channel?impulse?response?%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if?（Rayleigh）
????N=50;
????fm=100;
????B=20e3;
????fd=（rand（1N）-0.5）*2*fm;
????theta=randn（1N）*2*pi;
????c=randn（1N）;
????c=c/sum（c.^2）;
????t=0:fm/B:10000*fm/B;
????Tc=zeros（size（t））;
????
????Ts=zeros（size（t））;
????for?k=1:N
???????Tc=c（k）*cos（2*pi*fd（k）*t+theta（k））+Tc;
???????Ts=c（k）*sin（2*pi*fd（k）*t+theta（k））+Ts;
????end
????r=ones（Mt*Mr1）*（Tc.^2+Ts.^2）.^0.5;
????index=floor（rand（Mt*MrDS）*5000+1）;
end

MEE1=zeros（1length（SNR_dB））;
MEE2=zeros（1length（SNR_dB））;

for?snrl=1:length（SNR_dB）
????snrl
????estimation_error1=zeros（Mt*MrNsc）;
????estimation_error2=zeros（Mt*MrNsc）;
????R1=besselj（02*pi*fm*（Nsc+Ng）/B）;
????sigma2=10^（-SNR_dB（snrl）/10）;
????aa=（1-R1^2）/（1-R1^2+sigma2）;
????bb=sigma2*R1/（1-R1^2+sigma2）;

????for?iteration=1:iteration_max
????????%iteration????
????????if?AWGN==1
????????????h=ones（Mt*Mr1）;
????????else
????????????phi=rand*2*pi;
????????????h=r（index+iteration）*exp（j*phi）;
????????????%h=rand（Mt*MrDS）;
????????????h=h.*（ones（Mt*Mr1）*（exp（-0.5）.^[1:DS]））;
????????????h=h./（sqrt（sum（abs（h）.^22））*ones（1DS））;
????????end



????????CL=size（h2）;???????????????????????????????????????????????%?channel?length
????????data_time=zeros（MtNsc+Ng）;
????????data_qam=zeros（MtNsc）;
????????data_out=zeros（MrNsc）;
????????output=zeros（MrNsc）;

????????for?tx=1:Mt
????????????data_b=0*round（rand（4Nsc））;??????????????????????????????????%?data
????????????data_qam（tx:）=j*（2*（mod（data_b（1:）+data_b（2:）2）+2*data_b（1:））-3）+...
????????????2*（mod（data_b（3:）+data_b（4:）2）+2*data_b（3:））-3;
????????????for?loop=1:Mt?
????????????????data_qam（txpilots+loop-1）=（1+j）*（loop==tx）;??????????????%?pilots
????????????end
????????????data_time_temp=ifft（data_qam（tx:））;
????????????data_time（tx:）=[data_time_temp（end-Ng+1:end）?data_time_temp];
????????end
????
????????for?rx=1:Mr
????????????for?tx=1:Mt
????????????????output_temp=conv（data_time（tx:）h（（rx-1）*Mt+tx:））;
????????????????output（rx:）=output_temp（Ng+1:Ng+Nsc）+output（rx:）;
????????????end
????????????np=（sum（abs（output（rx:））.^2）/length（output（rx:）））*sigma2;
????????????noise=（randn（size（