function [BER] = MQAM_Simulation(longitud,M,SNR)
close all;
% Setting up parameters.
% Adding zeros to amount of bits to do it 2*M multiple.
resto = rem(longitud,2*M);
% Creation of random bits
bin = randint(longitud+(2*M-resto),1,2);
L=length(bin);
% Definition of bits with enough time to be observed. Setting up variables.
k=100;
f=10;
bit1=ones(1,k);
bit0=0*bit1;
Nivel=M;
symbol=ones(1,2*Nivel*k);
mbit=[];mx=[];my=[];
bin = fliplr(bin);
% Modulation in 2M bits groups
x=0;
y=0;
% We create groups of 2M bits to modulate
for n=0:2*Nivel:L-2*Nivel
bit=[];
xi=0;
yi=0;
% We assign an amplitude both 'x' and 'y' to the group
for m= 1:2:2*Nivel
if bin(n+m)==0
xi=xi+(2^((m-1)/2));
bit=[bit bit0];
else
xi=xi-(2^((m-1)/2));
bit=[bit bit1];
end
if bin(n+m+1)==0
yi=yi+(2^((m-1)/2));
bit=[bit bit0];
else
yi=yi-(2^((m-1)/2));
bit=[bit bit1];
end
end
x=xi*symbol;
y=yi*symbol;
% We store the generated symbol with the calculated amplitude inside mx and
% my variables. We update the mbit string with the last 2M modulated bits.
mx=[mx x];
my=[my y];
mbit=[mbit bit];
end
% Modulating and showing the constellation inside a scatterplot.
v=0:2*pi/k:2*pi*L-2*pi/k;
msync = mx + my*j;
qam=real(msync).*cos(f*v)-imag(msync).*sin(f*v);
scatterplot(msync),grid,xlabel('I'),ylabel('Q'),title('Constellation before sending');
pause;
% Addition of AWGB to simulate the channel
Vn=awgn(qam,SNR,'measured');
% Demodulation of received signal
Vnx=Vn.*cos(f*v);
Vny=Vn.*-1.*sin(f*v);
% Low pass filtering with a Butterworth filter
[b,a]=butter(2,0.04);
Hx=2.*filter(b,a,Vnx);
Hy=2.*filter(b,a,Vny);
ML=length(Hx);
% Showing the constellation received
msync2=[];
for m=k:2*Nivel*k:ML
Haux = Hx(m) + Hy(m)*j;
msync2 = [msync2 Haux];
end
scatterplot(msync2),grid,xlabel('I'),ylabel('Q'),title('Constellation received');
pause;
mdeb=[];
for m=k:2*Nivel*k:ML
sym=[];
thx=0;thy=0;
fx=0;fy=0;
for n=1:Nivel
if Hy(m) > thy
sym = [sym bit0];
fy=1;
else
sym = [sym bit1];
fy=-1;
end
if Hx(m) > thx
sym = [sym bit0];
fx=1;
else
sym = [sym bit1];
fx=-1;
end
thy = thy + fy*(2^(Nivel-n));
thx = thx + fx*(2^(Nivel-n));
end
mdeb=[mdeb fliplr(sym)];
end
% Calculation of wrong bits
err = 0;
bin2 = [];
for n=1:length(bin)
if(bin(n) == 1)
bin2 = [bin2 bit1];
elseif bin(n) == 0
bin2 = [bin2 bit0];
end
end
for n=1:length(bin2)
if(bin2(n) ~= mdeb(n))
err = err + 1;
end
end
error_bits_total = round(err/k);
% Calculation of BER to return the result
ber = error_bits_total/length(bin);
BER = ber;
% Showing final results
disp(['Total wrong bits = ' num2str(error_bits_total)]);
disp(['BER = ' num2str(ber)]);
figure(1);
subplot(4,1,1),plot(mbit,'r','linewidth',2),axis([0 k*L -0.5 1.5]),grid on,legend('Input data');
subplot(4,1,2),plot(qam,'m','linewidth',1.5),axis([0 k*L min(qam)*1.1 max(qam)*1.1]),grid on,legend('QAM Modulation');
subplot(4,1,3),plot(Vn,'g','linewidth',1.5),axis([0 k*L min(Vn)*1.1 max(Vn)*1.1]),grid on,legend('QAM Modulation with AWGN');
subplot(4,1,4),plot(mdeb,'k','linewidth',1.5),axis([0 k*L -0.5 1.5]),grid on,legend('Output data');
