%%*********************************************************************
%%%% 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=resha
