clear;clc;
% time vary channel
% 1 pilot block and 4 data block in the same coherent channel(frame=5)
% use 3 order polynomial h(t)=at^3+bt^2+ct+d
errorprob=zeros(1,6);
Nt=4;Nr=4;% number of Tx anntena and Rx anntena
% pilot matrix for 3bits BPSK 4by4 ,and normalize to Ep=1;
xp_bit=[0 0 1 1;0 1 0 1;0 0 0 0];% demap
symBitNum=1;
xp=eye(Nt);
Ep=trace(xp*xp')./Nt;
Es=1;
SNRScale=4;
frame=5;
blocksize=4;
order=3;
addpath('D:\ӤG\Papers\jakes_model');
eH=zeros(4,6); %row:SNR, col:veloctiy
for velocity=30:20:90;% km/hr
% velocity=50;
%%
for SNR_dB=20:SNRScale:20 % Eb/N0
% SNR_dB=50;
SNR=symBitNum*10^(SNR_dB/10);% Eb/N0
% forgetFactor=0.77^(SNR_dB/SNRScale+1+(round(velocity./10)-1)/2);
display(['v= ' num2str(velocity) ', SNR_dB= ' num2str(SNR_dB)]);
timeIter=0;
errorBitsCount=0;
totalBit=0;
totIter=0;%count the number of blocks
errorMax=10^7;
haha=0;
h=jakes_me(velocity,Nt,Nr,blocksize);
% total channel time =0.1s, sampling time =Nt*0.1ms, symbol time=0.1ms ,same H coef in a block
y=[];
x_bit_memo=[];%store trasnmitted signal bits
Hfinal=[];
while(errorBitsCount<errorMax)
if errorBitsCount>haha+499
haha=errorBitsCount;
display(['errorBitsCount= ' num2str(errorBitsCount) ' totIter= ' num2str(totIter)]);
end
sizey=size(y);
while (sizey(2)<4*(order*frame+1)) % at least tansmited (order-1) plitot blocks to solve the polynomial factor,which equals to order*frame+1 time slots
timeIter=timeIter+1;
% for the first time use pilot
if mod(timeIter,frame)==1
% channel generation
H=reshape(h(:,timeIter),Nr,Nt);
n=sqrt(Ep/(2*SNR))*randn(Nr,blocksize)+1j*sqrt(Ep/(2*SNR))*randn(Nr,blocksize);
y=[y H*xp+n]; % received pilot
timeIter=timeIter+1;
end
%check if channel length is enougth
if timeIter>length(h)
%% Count Channel MSE
long=length(Hfinal);%%%%%%%frame5 20,40 (240)
eHTmp=sum(sum(abs(Hfinal-h(:,1:long)))); %1 set ok
eH(velocity./20-0.5,SNR_dB/SNRScale+1)=eHTmp+eH(velocity./20-0.5,SNR_dB/SNRScale+1);
totIter=totIter+long;%%%%%%%%%%% frame5O2(-2),frame20O2,40O2,frame5O3(-12) ݳ̫_
if totIter==long*5 %check if there is 5 sets
eH(velocity./20-0.5,SNR_dB/SNRScale+1)=eH(velocity./20-0.5,SNR_dB/SNRScale+1)./(5)./(long*16);
errorBitsCount=errorMax+1;
Hfinal=[];
break;
end
Hfinal=[];
h=jakes_me(velocity,Nt,Nr,blocksize);
timeIter=0;
y=[];
sizey=size(y);
x_bit_memo=[];
continue; %check for frame 40
end
sizey=size(y);
if (sizey(2)>=4*(order*frame+1))%check if it is the 4rd pilot block,tpilotn뺡p
timeIter=timeIter-1;
break;
end
% generate signal bit
x_bit=randsrc(3,Nt,[0 1]);
x_bit_memo=[x_bit_memo x_bit];
x=zeros(Nt,Nt);
% signal symbol mapping
for index1=1:Nt
transAnt=2*x_bit(1,:)+x_bit(2,:)+1;
x(transAnt(index1),index1)=-2*x_bit(3,index1)+1;%signal symbol mapping
end
%thourgh channel
H=reshape(h(:,timeIter),Nr,Nt);
n=sqrt(Es/(2*SNR))*randn(Nr,Nt)+1j*sqrt(Es/(2*SNR))*randn(Nr,Nt);
y=[y H*x+n];
sizey=size(y);
end
if errorBitsCount>errorMax
continue;
end
%% Channel Estimaiton
Hp_hat=reshape([y(:,1:4) y(:,1+4*frame:4+4*frame) y(:,1+2*4*frame:4+2*4*frame) y(:,1+3*4*frame:4+3*4*frame)],16,order+1);
T=[(timeIter-3*frame)^3 (timeIter-2*frame)^3 (timeIter-frame)^3 timeIter^3;
(timeIter-3*frame)^2 (timeIter-2*frame)^2 (timeIter-frame)^2 timeIter^2;
timeIter-3*frame timeIter-2*frame timeIter-frame timeIter;ones(1,order+1)];
t=(timeIter-order*frame):timeIter;%time last
C=(Hp_hat*inv(T)).';
H_hat=C.'*[t.^3;t.^2;t;ones(1,order*frame+1)];%estimated channel for a batch
Hfinal=[ Hfinal H_hat(:,1:order*frame)];%̫@pilot block s
%% Signal Detection
for iter=2:frame
[x_hat x_hat_bit]=ML_SM_4x4_BPSK(y(:,4*iter-3:4*iter),reshape(H_hat(:,iter),4,4));
errorBitsCount=sum(sum(x_hat_bit~=x_bit_memo(:,4*(iter-1)-3:4*(iter-1))))+errorBitsCount;
end
for iter2=frame+2:2*frame
[x_hat x_hat_bit]=ML_SM_4x4_BPSK(y(:,4*iter2-3:4*iter2),reshape(H_hat(:,iter2),4,4));
errorBitsCount=sum(sum(x_hat_bit~=x_bit_memo(:,4*(iter2-2)-3:4*(iter2-2))))+errorBitsCount;
end
for iter3=(order-1)*frame+2:order*frame
[x_hat x_hat_bit]=ML_SM_4x4_BPSK(y(:,4*iter3-3:4*iter3),reshape(H_hat(:,iter3),4,4));
errorBitsCount=sum(sum(x_hat_bit~=x_bit_memo(:,4*(iter3-order)-3:4*(iter3-order))))+errorBitsCount;
end
%% empty the y and x
ytemp=y;
y=[];
y=ytemp(:,end-3:end);
x_bit_memo=[];
end
[row,col] = size(x_bit);
totalBit=row.*col*(totIter+timeIter-ceil(totIter./frame)-ceil(timeIter./frame));%totDataBit
errorprob(SNR_dB/SNRScale+1)=errorBitsCount./totalBit;
end
errorprobVel(velocity./20-0.5,:)=errorprob;
end
% SNR_dB=0:SNRScale:20;
% semilogy(SNR_dB,errorprobVel.','^-');
% grid on;
% legend('30 Km/hr','50 Km/hr','70 Km/hr','90 Km/hr','Location','best');
% % legend('70 Km/hr','90 Km/hr','Location','best');
% xlabel('Eb/No');ylabel('BER');
SNR_dB=0:SNRScale:20;
semilogy(SNR_dB,eH.','^-');
grid on;
legend('30 Km/hr','50 Km/hr','70 Km/hr','90 Km/hr','Location','best');
xlabel('Eb/No');ylabel('MSE');
% subplot(2,1,1)
% plot(real(Hfinal(10,:)),'o-r');
% hold on;
% plot(real(h(10,:)),'o-');
% grid on;
% legend('Estimated Channel','Reality Channel','Location','best');
% xlabel('Time(x0.4ms)');
% title(['Velocity ' num2str(velocity) ' Km/hr Eb/No=' num2str(SNR_dB) ' dB Real part']);
%
% subplot(2,1,2)
% plot(imag(Hfinal(10,:)),'o-r')
% hold on;
% plot(imag(h(10,:)),'o-')
% grid on;
% legend('Estimated Channel','Reality Channel','Location','best');
% title(['Velocity ' num2str(velocity) ' Km/hr Eb/No=' num2str(SNR_dB) ' dB Imag part']);
% xlabel('Time(x0.4ms)')
%
