clear,clf
N = 30000;
t1=cputime;
Rb = 10000;
Tb = 1/Rb; % bit symbol
fc = 4*Rb; % tan so song mang
fs = 16*Rb; % tan so lay mau
Ts = 1/fs;
M1 = 2; k1 = log2(M1);
M2 = 4; k2 = log2(M2); L2 = 2^(k2/2);
M3 = 16; k3 = log2(M3); L3 = 2^(k3/2);
M4 = 64; k4 = log2(M4); L4 = 2^(k4/2);
% kenh truyen fading
v = 40; % v = 20~80km/h
fd = 1000*v*fc/(3e8);
chan1 = rayleighchan(1/Rb,fd); %100);
chan1.PathDelays = [0, 1.5*1e-6, 2.3*1e-6, 2.69*1e-6];
chan1.AvgPathGaindB = [-1.4, -8.08, -11.42, -13.09];
chan1.NormalizePathGains = 1;
chan1.ResetBeforeFiltering = 0;
% chan1.StoreHistory = true;
% kenh truyen hilly terrain
chan2 = rayleighchan(1/Rb,fd);% 100);
chan2.PathDelays = [0, 0.9*1e-6, 1.8*1e-6, 15.04*1e-6];
chan2.AvgPathGaindB = [-0.7, -14.99, -15.44, -29.30];
chan2.NormalizePathGains = 1;
chan2.ResetBeforeFiltering = 0;
% khoi tao du lieu
data = randint(1,N);
% dieu che bpsk
tx_bpsk = codebpsk(data,M1,Rb,fc,fs);
fadbpsk1 = abs(filter(chan1,ones(size(tx_bpsk))));
fadbpsk2 = abs(filter(chan2,ones(size(tx_bpsk))));
bpskFad1 = fadbpsk1.*tx_bpsk;
bpskFad2 = fadbpsk2.*tx_bpsk;
% dieu che qpsk
tx_qpsk = codeqpsk(data,M2,Rb,fc,fs);
fadqpsk1 = abs(filter(chan1,ones(size(tx_qpsk))));
fadqpsk2 = abs(filter(chan2,ones(size(tx_qpsk))));
qpskFad1 = fadqpsk1.*tx_qpsk;
qpskFad2 = fadqpsk2.*tx_qpsk;
% dieu che 16qam
tx_qam16 = codeqam(data,M3,Rb,fc,fs);
fadqam161 = abs(filter(chan1,ones(size(tx_qam16))));
fadqam162 = abs(filter(chan2,ones(size(tx_qam16))));
qam16Fad1 = fadqam161.*tx_qam16;
qam16Fad2 = fadqam162.*tx_qam16;
% dieu che 64qam
tx_qam64 = codeqam(data,M4,Rb,fc,fs);
fadqam641 = abs(filter(chan1,ones(size(tx_qam64))));
fadqam642 = abs(filter(chan2,ones(size(tx_qam64))));
qam64Fad1 = fadqam641.*tx_qam64;
qam64Fad2 = fadqam642.*tx_qam64;
% ve tin hieu
% subplot(4,1,1); plot(tx_bpsk); title('dieu che bpsk'); grid on; axis on;
% subplot(4,1,2); plot(tx_qpsk); title('dieu che qpsk'); grid on; axis on;
% subplot(4,1,3); plot(tx_qam16); title('dieu che 16qam'); grid on; axis on;
% subplot(4,1,4); plot(tx_qam64); title('dieu che 64qam'); grid on; axis on;
EbN0dB = 1:30; % energy of bit per noise (dB)
EbN0 = 10.^(EbN0dB/10);
% vong lap theo tung gia tri EbN0dB
for i=1:length(EbN0dB)
for inter = 1:5 % so lan lap cho mot gia tri EbN0
% nhieu awgn cho dieu che bpsk
Varbpsk = sqrt(1./(2*(log2(M1))*EbN0(i)));
noisebpsk = Varbpsk*(randn(1,length(bpskFad1)) + 1i*randn(1,length(bpskFad1)));
t_bpsk1 = bpskFad1 + noisebpsk; %awgn(bpskFad1,EbN0dB(i),'measured');%
t_bpsk2 = bpskFad2 + noisebpsk; %awgn(bpskFad2,EbN0dB(i),'measured');% + noisebpsk;
t_bpsk = tx_bpsk + noisebpsk; %awgn(tx_bpsk,EbN0dB(i),'measured');% + noisebpsk;
rx_bpsk1 = codebpskdemod((t_bpsk1./fadbpsk1),fc,fs,Rb,N);
rx_bpsk2 = codebpskdemod((t_bpsk2./fadbpsk2),fc,fs,Rb,N);
rx_bpsk = codebpskdemod(t_bpsk,fc,fs,Rb,N);
% tinh ber
nebpsk1(i) = biterr(data,rx_bpsk1);
nebpsk2(i) = biterr(data,rx_bpsk2);
nebpsk(i) = biterr(data,rx_bpsk);
berbpsk(i) = nebpsk(i)/N;
berbpsk1(i) = nebpsk1(i)/N;
berbpsk2(i) = nebpsk2(i)/N;
% nhieu awgn cho dieu che qpsk
Varqpsk = sqrt(1./(2*(log2(M2))*EbN0(i)));
noiseqpsk = Varqpsk*(randn(1,length(qpskFad1)) + 1i*randn(1,length(qpskFad1)));
t_qpsk = tx_qpsk + noiseqpsk; %awgn(tx_qpsk,EbN0dB(i),'measured');%
t_qpsk1 = qpskFad1 + noiseqpsk; % awgn(qpskFad1,EbN0dB(i),'measured');%
t_qpsk2 = qpskFad2 + noiseqpsk; % awgn(qpskFad2,EbN0dB(i),'measured');%
rx_qpsk = codeqpskdemod(t_qpsk,fc,fs,Rb,N);
rx_qpsk1 = codeqpskdemod((t_qpsk1./fadqpsk1),fc,fs,Rb,N);
rx_qpsk2 = codeqpskdemod((t_qpsk2./fadqpsk2),fc,fs,Rb,N);
% tinh ber
neqpsk(i) = biterr(data,rx_qpsk);
neqpsk1(i) = biterr(data,rx_qpsk1);
neqpsk2(i) = biterr(data,rx_qpsk2);
berqpsk(i) = neqpsk(i)/N;
berqpsk1(i) = neqpsk1(i)/N;
berqpsk2(i) = neqpsk2(i)/N;
% nhieu awgn dieu che qam16
varqam16 = sqrt(1./(2*(log2(M3))*EbN0(i)));
noiseqam16 = varqam16*(randn(1,length(qam16Fad1)) + 1i*randn(1,length(qam16Fad1)));
t_qam16 = tx_qam16 + noiseqam16; % awgn(tx_qam16,EbN0dB(i),'measured');%
t_qam161 = qam16Fad1 + noiseqam16; % awgn(qam16Fad1,EbN0dB(i),'measured'); %
t_qam162 = qam16Fad2 + noiseqam16; % awgn(qam16Fad2,EbN0dB(i),'measured'); %
rx_qam16 = codeqamde(t_qam16,fc,fs,Rb,N,M3);
rx_qam161 = codeqamde((t_qam161./fadqam161),fc,fs,Rb,N,M3);
rx_qam162 = codeqamde((t_qam162./fadqam162),fc,fs,Rb,N,M3);
% tinh ber
neqam16(i) = biterr(data,rx_qam16);
neqam161(i) = biterr(data,rx_qam161);
neqam162(i) = biterr(data,rx_qam162);
berqam16(i) = neqam16(i)/N;
berqam161(i) = neqam161(i)/N;
berqam162(i) = neqam162(i)/N;
% nhieu awgn dieu che qam64
varqam64 = sqrt(1./(2*(log2(M4))*EbN0(i)));
noiseqam64 = varqam64*(randn(1,length(qam64Fad1)) + 1i*randn(1,length(qam64Fad1)));
t_qam64 = tx_qam64 + noiseqam64;
t_qam641 = qam64Fad1 + noiseqam64;
t_qam642 = qam64Fad2 + noiseqam64;
rx_qam64 = codeqamde(t_qam64,fc,fs,Rb,N,M4);
rx_qam641 = codeqamde((t_qam641./fadqam641),fc,fs,Rb,N,M4);
rx_qam642 = codeqamde((t_qam642./fadqam642),fc,fs,Rb,N,M4);
% tinh ber
neqam64(i) = biterr(data,rx_qam64);
neqam641(i) = biterr(data,rx_qam641);
neqam642(i) = biterr(data,rx_qam642);
berqam64(i) = neqam64(i)/N;
berqam641(i) = neqam641(i)/N;
berqam642(i) = neqam642(i)/N;
end
end
% tinh ber theo ly thuyet
bpskFadTh = berfading(EbN0dB,'psk',M1,1);
qpskFadTh = berfading(EbN0dB,'psk',M2,1);
qam16FadTh = berfading(EbN0dB,'qam',M3,1);
qam64FadTh = berfading(EbN0dB,'qam',M4,1);
% bpskTh = berawgn(EbN0dB,'psk',M1);
% ve gian do ber - ebn0
% figure;
semilogy(EbN0dB,berbpsk,'b-',EbN0dB,berbpsk1,'b-o',EbN0dB,berbpsk2,'b-x',...
EbN0dB,bpskFadTh,'r-o',EbN0dB,berqpsk,'g-',EbN0dB,berqpsk1,'g-o',EbN0dB,...
berqpsk2,'g-x',EbN0dB,qpskFadTh,'g-.',EbN0dB,berqam16,'k-',EbN0dB,berqam161,'k-o',EbN0dB,...
berqam162,'k-x',EbN0dB,qam16FadTh,'k-.',EbN0dB,berqam64,'m-',EbN0dB,berqam641,'m-o',EbN0dB,...
berqam642,'m-x',EbN0dB,qam64FadTh,'m-.' );
legend('bpsk awgn','bpsk chan1','bpsk chan2','bpsk fading',....
'qpsk awgn','qpsk chan1','qpsk chan2','qpsk theory fad',....
'qam16 awgn','qam16 chan1','qam16 chan2',...
'qam64 awgn','qam64 chan1','qam64 chan2');
axis([0 30 10^-5 1]);
grid on;
t2=cputime;
ttt = t2-t1;