%%*********************************************************************
%%%%     M.1225 信道下，MIMO,2transmitters and 2receivers ,進行ml估計、
%%%********************************************************************
function [Dem1,Dem2]=twoto22(out1,out2,SNR)
%=================
% common settings
%=================
wordsize=2;               % represent the mode of modulation, wordsize=2:QPSK, wordsize=4:16QAM, wordsize=6:64QAM;     
NumCarr=256;              % Number of transmission carriers    
Lcp=64;                   % Guard Interval(length of cyclic extension)=12.5% of NumCarr
FrameGuard=NumCarr+Lcp/2;   % Guard Time between successive frames (one symbol period)
Numsymb=100;              % the number of symbols
seed=1234;
rand('seed',seed);	      % Set to new seed
Wk=[1 0 1 0 0 0 1 1 1 0 1 0 0 0 1 1 1 0 1 0 0 0 1 1];
Pilot1_number=24;
Pilot2_number=24;
PilotIndex1=[37:8:221];     % pilot interval=24
PilotIndex2=[40:8:224];

fc=3.5e9;                 % operation frequency
V=75;                     % moving speed, V km/h
fdmax=V*fc/3e8/3.6;        % maximum frequency shift: fmax=V*fc/C

%generate s0
%Generate the pattern of data and pilots
Pattern=ones(NumCarr,Numsymb);    % the position of data is set as 1.
Pattern([1:28,230:256],:)=2;      % guard band=0
Pattern([37:8:221],:)=4;         % the position of pilot1 is set as 4.
Pattern([40:8:224],:)=5;          % the position of pilot5 is set as 5.
Pattern(129,:)=3;                 % DC=0

% Generate data and pilots
Data0=zeros(size(Pattern));
for n=1:Numsymb
    Data0([37:8:221],n)=2*(1-2*Wk).';       
end;
     
Data_Pattern=find(Pattern==1);                    %the pattern of data
Baohu_Pattern=find(Pattern==2);
Dc_Pattern=find(Pattern==3);
Pilot1_Pattern=find(Pattern==4);
Pilot2_Pattern=find(Pattern==5);
clear Pattern

Datatx0=floor(rand(1,length(Data_Pattern))*(2^wordsize));
bbb=zeros(1,500*30);
 for i=1:500*30
      aaa=num2str(out1(2*i-1:i*2)');
      bbb(i)=bin2dec(aaa');
  end;
Datatx0(1:500*30)=bbb;% Generate the data
Baohutx0=floor(rand(1,length(Baohu_Pattern))*(2^wordsize));

Tx0=dec2bin(Datatx0,wordsize)-48;
% Mapping to the signal constellation follow
mapping=get80216map(2^wordsize);% 調用.p文件。
for k=1:length(Datatx0)
    ModSignal0(k)=mapping(1+Datatx0(k));
end;
Data0(Data_Pattern)=ModSignal0;

for k=1:length(Baohutx0)
    BaohuSignal0(k)=mapping(1+Baohutx0(k));
end;
Data0(Baohu_Pattern)=BaohuSignal0;
Data0(Dc_Pattern)=1;
clear ModSignal0;
clear BaohuSigna0;

%===================================
% Find the time waveform using IFFT
%===================================
BaseSig0=ifft(Data0);      % Generating baseband signal s0 using IFFT 
BaseSig3=ifft(conj(Data0));

%generate s1
%Generate the pattern of data and pilots
Pattern=ones(NumCarr,Numsymb);    % the position of data is set as 1.
Pattern([1:28,230:256],:)=2;      % guard band=0
Pattern([37:8:221],:)=4;         % the position of pilot1 is set as 4.
Pattern([40:8:224],:)=5;          % the position of pilot5 is set as 5.
Pattern(129,:)=3;                 % DC=0

% Generate data and pilots
Data1=zeros(size(Pattern));
for n=1:Numsymb
    Data1([40:8:224],n)=2*(1-2*Wk).'; 
end;

Data_Pattern=find(Pattern==1);                    %the pattern of data
Baohu_Pattern=find(Pattern==2);
Dc_Pattern=find(Pattern==3);
Pilot1_Pattern=find(Pattern==4);
Pilot2_Pattern=find(Pattern==5);
clear Pattern

Datatx1=floor(rand(1,length(Data_Pattern))*(2^wordsize)); 
bbb=zeros(1,500*30);
 for i=1:500*30
      aaa=num2str(out2(2*i-1:i*2)');
      bbb(i)=bin2dec(aaa');
  end;
Datatx1(1:500*30)=bbb;% Generate the data
Baohutx1=floor(rand(1,length(Baohu_Pattern))*(2^wordsize));

Tx1=dec2bin(Datatx1,wordsize)-48;
% Mapping to the signal constellation follow
mapping=get80216map(2^wordsize);% 調用.p文件。
for k=1:length(Datatx1)
    ModSignal1(k)=mapping(1+Datatx1(k));
end;
Data1(Data_Pattern)=ModSignal1;

for k=1:length(Baohutx1)
    BaohuSignal1(k)=mapping(1+Baohutx1(k));
end;
Data1(Baohu_Pattern)=BaohuSignal1;
Data1(Dc_Pattern)=1;
clear ModSignal
clear BaohuSignal

%===================================
% Find the time waveform using IFFT
%===================================
BaseSig1=ifft(Data1);      % Generating baseband signal s1 using IFFT 
BaseSig2=ifft(-conj(Data1));


%=======================
% Adding Guard Interval
%=======================
BaseSignal0=[BaseSig0((end-(Lcp/2)+1):end,:); BaseSig0];
BaseSignal1=[BaseSig1((end-(Lcp/2)+1):end,:); BaseSig1];
BaseSignal2=[BaseSig2((end-(Lcp/2)+1):end,:); BaseSig2];
BaseSignal3=[BaseSig3((end-(Lcp/2)+1):end,:); BaseSig3];
	
% =======================
% M.1225 channel 
% =======================
fade=Rayleigh(fdmax);
%1st transmitter to 1st receiver H11*BaseSignal0
path11_1=ones(FrameGuard,1)*fade(1,[1:Numsymb]).*BaseSignal0;
path11_2=ones(FrameGuard,1)*fade(2,[1:Numsymb]).*BaseSignal0;
path11_3=ones(FrameGuard,1)*fade(3,[1:Numsymb]).*BaseSignal0;
path11_4=ones(FrameGuard,1)*fade(4,[1:Numsymb]).*BaseSignal0;
path11_5=ones(FrameGuard,1)*fade(5,[1:Numsymb]).*BaseSignal0;
path11_6=ones(FrameGuard,1)*fade(6,[1:Numsymb]).*BaseSignal0;
path11_1=reshape(path11_1,1,size(path11_1,1)*size(path11_1,2));
path11_2=reshape(path11_2,1,size(path11_2,1)*size(path11_2,2));
path11_3=reshape(path11_3,1,size(path11_3,1)*size(path11_3,2));
path11_4=reshape(path11_4,1,size(path11_4,1)*size(path11_4,2));
path11_5=reshape(path11_5,1,size(path11_5,1)*size(path11_5,2));
path11_6=reshape(path11_6,1,size(path11_6,1)*size(path11_6,2));
%the delay is [0 0.31 0.71 1.09 1.73 2.51] us
path11_1=[path11_1 zeros(1,10)];    % the largest delay is 10 sample                                                                                                                          
path11_2=[zeros(1,1) path11_2 zeros(1,9)];
path11_3=[zeros(1,3) path11_3 zeros(1,7)];
path11_4=[zeros(1,4) path11_4 zeros(1,6)];
path11_5=[zeros(1,7) path11_5 zeros(1,3)];
path11_6=[zeros(1,10) path11_6];
RxSignal00=path11_1+path11_2+path11_3+path11_4+path11_5+path11_6;
RxSignal00=RxSignal00(1:FrameGuard*Numsymb);

%2st transmitter to 1st receiver H21*BaseSignal1
path21_1=ones(FrameGuard,1)*fade(1,[1+1000:Numsymb+1000]).*BaseSignal1;
path21_2=ones(FrameGuard,1)*fade(2,[1+1000:Numsymb+1000]).*BaseSignal1;
path21_3=ones(FrameGuard,1)*fade(3,[1+1000:Numsymb+1000]).*BaseSignal1;
path21_4=ones(FrameGuard,1)*fade(4,[1+1000:Numsymb+1000]).*BaseSignal1;
path21_5=ones(FrameGuard,1)*fade(5,[1+1000:Numsymb+1000]).*BaseSignal1;
path21_6=ones(FrameGuard,1)*fade(6,[1+1000:Numsymb+1000]).*BaseSignal1;
path21_1=reshape(path21_1,1,size(path21_1,1)*size(path21_1,2));
path21_2=reshape(path21_2,1,size(path21_2,1)*size(path21_2,2));
path21_3=reshape(path21_3,1,size(path21_3,1)*size(path21_3,2));
path21_4=reshape(path21_4,1,size(path21_4,1)*size(path21_4,2));
path21_5=reshape(path21_5,1,size(path21_5,1)*size(path21_5,2));
path21_6=reshape(path21_6,1,size(path21_6,1)*size(path21_6,2));
path21_1=[path21_1 zeros(1,10)];    % the largest delay is 10 sample                                                                                                                          
path21_2=[zeros(1,1) path21_2 zeros(1,9)];
path21_3=[zeros(1,3) path21_3 zeros(1,7)];
path21_4=[zeros(1,4) path21_4 zeros(1,6)];
path21_5=[zeros(1,7) path21_5 zeros(1,3)];
path21_6=[zeros(1,10) path21_6];
RxSignal11=path21_1+path21_2+path21_3+path21_4+path21_5+path21_6;
RxSignal11=RxSignal11(1:FrameGuard*Numsymb);

%1st transmitter to 1st receiver H11*BaseSignal2
path11_1=ones(FrameGuard,1)*fade(1,[1:Numsymb]).*BaseSignal2;
path11_2=ones(FrameGuard,1)*fade(2,[1:Numsymb]).*BaseSignal2;
path11_3=ones(FrameGuard,1)*fade(3,[1:Numsymb]).*BaseSignal2;
path11_4=ones(FrameGuard,1)*fade(4,[1:Numsymb]).*BaseSignal2;
path11_5=ones(FrameGuard,1)*fade(5,[1:Numsymb]).*BaseSignal2;
path11_6=ones(FrameGuard,1)*fade(6,[1:Numsymb]).*BaseSignal2;
path11_1=reshape(path11_1,1,size(path11_1,1)*size(path11_1,2));
path11_2=reshape(path11_2,1,size(path11_2,1)*size(path11_2,2));
path11_3=reshape(path11_3,1,size(path11_3,1)*size(path11_3,2));
path11_4=reshape(path11_4,1,size(path11_4,1)*size(path11_4,2));
path11_5=reshape(path11_5,1,size(path11_5,1)*size(path11_5,2));
path11_6=reshape(path11_6,1,size(path11_6,1)*size(path11_6,2));
path11_1=[path11_1 zeros(1,10)];    % the largest delay is 10 sample                                                                                                                          
path11_2=[zeros(1,1) path11_2 zeros(1,9)];
path11_3=[zeros(1,3) path11_3 zeros(1,7)];
path11_4=[zeros(1,4) path11_4 zeros(1,6)];
path11_5=[zeros(1,7) path11_5 zeros(1,3)];
path11_6=[zeros(1,10) path11_6];
RxSignal01=path11_1+path11_2+path11_3+path11_4+path11_5+path11_6;
RxSignal01=RxSignal01(1:FrameGuard*Numsymb);

%2st transmitter to 1st receiver H21*BaseSignal3
path21_1=ones(FrameGuard,1)*fade(1,[1+1000:Numsymb+1000]).*BaseSignal3;
path21_2=ones(FrameGuard,1)*fade(2,[1+1000:Numsymb+1000]).*BaseSignal3;
path21_3=ones(FrameGuard,1)*fade(3,[1+1000:Numsymb+1000]).*BaseSignal3;
path21_4=ones(FrameGuard,1)*fade(4,[1+1000:Numsymb+1000]).*BaseSignal3;
path21_5=ones(FrameGuard,1)*fade(5,[1+1000:Numsymb+1000]).*BaseSignal3;
path21_6=ones(FrameGuard,1)*fade(6,[1+1000:Numsymb+1000]).*BaseSignal3;
path21_1=reshape(path21_1,1,size(path21_1,1)*size(path21_1,2));
path21_2=reshape(path21_2,1,size(path21_2,1)*size(path21_2,2));
path21_3=reshape(path21_3,1,size(path21_3,1)*size(path21_3,2));
path21_4=reshape(path21_4,1,size(path21_4,1)*size(path21_4,2));
path21_5=reshape(path21_5,1,size(path21_5,1)*size(path21_5,2));
path21_6=reshape(path21_6,1,size(path21_6,1)*size(path21_6,2));
path21_1=[path21_1 zeros(1,10)];    % the largest delay is 10 sample                                                                                                                          
path21_2=[zeros(1,1) path21_2 zeros(1,9)];
path21_3=[zeros(1,3) path21_3 zeros(1,7)];
path21_4=[zeros(1,4) path21_4 zeros(1,6)];
path21_5=[zeros(1,7) path21_5 zeros(1,3)];
path21_6=[zeros(1,10) path21_6];
RxSignal10=path21_1+path21_2+path21_3+path21_4+path21_5+path21_6;
RxSignal10=RxSignal10(1:FrameGuard*Numsymb);

%1st transmitter to 2st receiver H12*BaseSignal0
path12_1=ones(FrameGuard,1)*fade(1,[1+2000:Numsymb+2000]).*BaseSignal0;
path12_2=ones(FrameGuard,1)*fade(2,[1+2000:Numsymb+2000]).*BaseSignal0;
path12_3=ones(FrameGuard,1)*fade(3,[1+2000:Numsymb+2000]).*BaseSignal0;
path12_4=ones(FrameGuard,1)*fade(4,[1+2000:Numsymb+2000]).*BaseSignal0;
path12_5=ones(FrameGuard,1)*fade(5,[1+2000:Numsymb+2000]).*BaseSignal0;
path12_6=ones(FrameGuard,1)*fade(6,[1+2000:Numsymb+2000]).*BaseSignal0;
path12_1=reshape(path12_1,1,size(path12_1,1)*size(path12_1,2));
path12_2=reshape(path12_2,1,size(path12_2,1)*size(path12_2,2));
path12_3=reshape(path12_3,1,size(path12_3,1)*size(path12_3,2));
path12_4=reshape(path12_4,1,size(path12_4,1)*size(path12_4,2));
path12_5=reshape(path12_5,1,size(path12_5,1)*size(path12_5,2));
path12_6=reshape(path12_6,1,size(path12_6,1)*size(path12_6,2));
path12_1=[path12_1 zeros(1,10)];    % the largest delay is 10 sample                                                                                                                          
path12_2=[zeros(1,1) path12_2 zeros(1,9)];
path12_3=[zeros(1,3) path12_3 zeros(1,7)];
path12_4=[zeros(1,4) path12_4 zeros(1,6)];
path12_5=[zeros(1,7) path12_5 zeros(1,3)];
path12_6=[zeros(1,10) path12_6];
RxSignal20=path12_1+path12_2+path12_3+path12_4+path12_5+path12_6;
RxSignal20=RxSignal20(1:FrameGuard*Numsymb);

%2st transmitter to 2st receiver H22*BaseSignal1
path22_1=ones(FrameGuard,1)*fade(1,[1+3000:Numsymb+3000]).*BaseSignal1;
path22_2=ones(FrameGuard,1)*fade(2,[1+3000:Numsymb+3000]).*BaseSignal1;
path22_3=ones(FrameGuard,1)*fade(3,[1+3000:Numsymb+3000]).*BaseSignal1;
path22_4=ones(FrameGuard,1)*fade(4,[1+3000:Numsymb+3000]).*BaseSignal1;
path22_5=ones(FrameGuard,1)*fade(5,[1+3000:Numsymb+3000]).*BaseSignal1;
path22_6=ones(FrameGuard,1)*fade(6,[1+3000:Numsymb+3000]).*BaseSignal1;
path22_1=reshape(path22_1,1,size(path22_1,1)*size(path22_1,2));
path22_2=reshape(path22_2,1,size(path22_2,1)*size(path22_2,2));
path22_3=reshape(path22_3,1,size(path22_3,1)*size(path22_3,2));
path22_4=reshape(path22_4,1,size(path22_4,1)*size(path22_4,2));
path22_5=reshape(path22_5,1,size(path22_5,1)*size(path22_5,2));
path22_6=reshape(path22_6,1,size(path22_6,1)*size(path22_6,2));
path22_1=[path22_1 zeros(1,10)];    % the largest delay is 10 sample                                                                                                                          
path22_2=[zeros(1,1) path22_2 zeros(1,9)];
path22_3=[zeros(1,3) path22_3 zeros(1,7)];
path22_4=[zeros(1,4) path22_4 zeros(1,6)];