Monte-Carlo Simulation of QPSK Communication System in MATLAB #Communication# Signal Processing #
I really did much programming and simulations in R, Python and C++ last semseter. Suddently, I found that I haven't been using MATLAB for almost two years! Though it is an easy language (assume it is a language, anyway) compared with...er like cpp, I'd better review some of the MATLAB project I did before.
Here, I am going to present a MATLAB based simulation of a QPSK communication system model with various source coding methods (none/gray/hamming code) and with different judegement method (min Euclidean/max projection). the results are followed
Main function:
****************************************************************
clc; clear; A = 1;L = 5000; o2 = 0.05:0.05:1; for j = 1:20 [ an,c1 ] = ty2_quanternary_isource( L ); [ an2 ] = ty2_4PSK_projection( an,A ); [ nt ] = ty2_gauss_random_number_generator( an,o2(j)); [ an3 ] = ty2_output_AWGN( an2,nt ); Es = 1; r = Es/(2*o2(j)); SNR1(j) = 10*log10(r); [ an4,c2 ] = ty2_min_Euclidean_distance_judgement(an3,an,A ); [ Pb1(j) ] = ty2_bit_error( c1,c2,an ); end for j = 1:20 [ an,c1 ] = ty2_quanternary_isource( L ); [ hm1] = ty2_hamming_code_isource( c1,L ); [ hm ] = ty2_hamming_hm( hm1 ); [ hm2 ] = ty2_4PSK_projection( hm,A ); [ nt ] = ty2_gauss_random_number_generator( hm,o2(j) ); [ hm3 ] = ty2_output_AWGN( hm2,nt ); Es = 1; r = Es/(2*o2(j)); SNR2(j) = 10*log10(r); [ n,hm4 ] = ty2_min_Euclidean_distance_judgement( hm3,hm,A ); [ hm5,hm6 ] = ty2_hamming_decode( hm4 ); [ c2 ] = ty2_Gray_code( hm6 ); [ Pb2(j) ] = ty2_bit_error( c1,c2,an ); end figure;semilogy(SNR1,Pb1,'g');hold on;grid on;semilogy(SNR2,Pb2,'r');
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~···
clc; clear; A = 1;L = 1000;o2 = [0,0.1,0.5,1]; for i = 1:4 [ an,c1 ] = ty2_quanternary_isource( L ); [ an2 ] = ty2_4PSK_projection( an,A ); [ nt ] = ty2_gauss_random_number_generator( an,o2(i) ); [ an3 ] = ty2_output_AWGN( an2,nt ); [ an4,c2 ] = ty2_min_Euclidean_distance_judgement(an3,an,A ); [ Pe ] = ty2_symbol_error( an,an4 ); [ Pb ] = ty2_bit_error( c1,c2,an ); scatterplot(an3); end
L = 5000; o2 = 0.05:0.05:1; for j = 1:20 [ an,c1 ] = ty2_quanternary_isource( L ); [ an2 ] = ty2_4PSK_projection( an,A ); [ nt ] = ty2_gauss_random_number_generator( an,o2(j)); [ an3 ] = ty2_output_AWGN( an2,nt ); Es = 1; r = Es/(o2(j)*2); SNR(j) = 10*log10(r); Q(j) = 0.5*erfc(sqrt(r)); [ an4,c2 ] = ty2_min_Euclidean_distance_judgement(an3,an,A ); [ Pb(j) ] = ty2_bit_error( c1,c2,an ); end
figure;semilogy(SNR,Pb,'g');hold on;grid on;semilogy(SNR,Q,'r'); xlabel('SNR'),ylabel('Pb');
4PSK Projection function:
************************************
function [ an2 ] = ty2_4PSK_projection( an,A ) %4PSK projection; an is source sequence while an2 is the projected sequence
L = length(an); an2 = zeros(L,2); for t = 1:L if an(t) == 0 an2(t,:) = [A,0]; else if an(t) == 1 an2(t,:) = [0,A]; else if an(t) == 2 an2(t,:) = [-A,0]; else if an(t) == 3 an2(t,:) = [0,-A]; end end end end end
end
Gaussian random number generator function:
******************************************************
function [ nt ] = ty2_gauss_random_number_generator( an,o2 ) %UNTITLED4 Summary of this function goes here % Detailed explanation goes here L = length(an); for i = 1:L u=rand; z=sqrt(o2)*sqrt(2*log(1/(1-u))); u=rand; nc(i)=z*cos(2*pi*u); ns(i)=z*sin(2*pi*u); end nt = [nc',ns']; end
Grey Gode function:
***************************************************
function [ c1 ] = ty2_Gray_code( an ) %an to be source code; c is its gray code L = length(an); c1 = zeros(L,2); for i = 1:L if an(i)==0 c1(i,:) = [0,0]; else if an(i)==1 c1(i,:) = [0,1]; else if an(i)==2 c1(i,:) = [1,1]; else if an(i)==3 c1(i,:) = [1,0]; end end end end end
end
Hamming code function:
***********************************************
function [ hm1] = ty2_hamming_code_isource( c1,L ) % hamming source c1 is gray code，hm1 is hamming code hm1 = zeros(L/2,7); for i = 1:L/2 hm1(i,[4:7]) = [c1(2*i-1,:),c1(2*i,:)]; hm1(i,1) = xor(hm1(i,7),xor(hm1(i,4),hm1(i,5))); hm1(i,2) = xor(hm1(i,7),xor(hm1(i,4),hm1(i,6))); hm1(i,3) = xor(hm1(i,7),xor(hm1(i,6),hm1(i,5))); end end
Hamming decoding function:
***************************************************************
function [ hm5,hm6 ] = ty2_hamming_decode( hm4 ) %hm4 is the input hamming code，hm5 is hamming code after error correction，hm6 is the correcting sequence; t = hm4'; t = t(:); L = length(t); hm5=zeros(L/7,7); for i = 1:L/7 hm5(i,:) = t(7*i-6:7*i); end for i = 1:L/7 s1 = xor( hm5(i,3),xor( hm5(i,5),xor( hm5(i,7), hm5(i,6)))); s2 = xor( hm5(i,2),xor( hm5(i,4),xor( hm5(i,7), hm5(i,6)))); s3 = xor( hm5(i,1),xor( hm5(i,4),xor( hm5(i,7), hm5(i,5)))); if [s1,s2,s3] == [0,1,1] hm5(i,4) = ~hm5(i,4); else if [s1,s2,s3] == [1,0,1] hm5(i,5) = ~hm5(i,5); else if [s1,s2,s3] == [1,1,0] hm5(i,6) = ~hm5(i,6); else if [s1,s2,s3] == [1,1,1] hm5(i,7) = ~hm5(i,7); end end end end end for i = 1:L/7 if hm5(i,4:5) ==[0,0] hm6(2*i-1) = 0; else if hm5(i,4:5) ==[0,1] hm6(2*i-1) = 1; else if hm5(i,4:5) ==[1,1] hm6(2*i-1) = 2; else if hm5(i,4:5) ==[1,0] hm6(2*i-1) = 3; end end end end if hm5(i,6:7) ==[0,0] hm6(2*i) = 0; else if hm5(i,6:7) ==[0,1] hm6(2*i) = 1; else if hm5(i,6:7) ==[1,1] hm6(2*i) = 2; else if hm5(i,6:7) ==[1,0] hm6(2*i) = 3; end end end end end
end
Max projection judegement function:
******************************************************************
function [ an4,c2 ] = ty2_max_projection_judgement( an3,an,A ) %，an3 is the received signal，an is source code，an4 is judged signal，c2 is current gray code L = length(an); for i = 1:L a(1) = an3(i,:)*[A,0]'; a(2) = an3(i,:)*[0,A]'; a(3) = an3(i,:)*[-A,0]'; a(4) = an3(i,:)*[0,-A]'; t = find(a==max(a)); an4(i) = t-1; end [ c2 ] = ty2_Gray_code( an4 );
end
Minimun Euclidean distance judgement function:
**********************************************************************
function [ an4,c2 ] = ty2_min_Euclidean_distance_judgement(an3,an,A ) L = length(an); for i = 1:L a(1) = (an3(i,:)-[A,0])*(an3(i,:)-[A,0])'; a(2) = (an3(i,:)-[0,A])*(an3(i,:)-[0,A])'; a(3) = (an3(i,:)-[-A,0])*(an3(i,:)-[-A,0])'; a(4) = (an3(i,:)-[0,-A])*(an3(i,:)-[0,-A])'; t = find(a==min(a)); an4(i) = t-1; end [ c2 ] = ty2_Gray_code( an4 );
end
