%  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                              
%                         CNTI
%   Communications Numériques et Théorie de l'Information
%
%     Fichier .m pour engendrer des diagrammes de l'oeil
%
%                     GRG  - 2009
%
%  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%     

% Nettoyage général avant de commencer
clear;
clf reset;

%  Variables de simulation
Nech_symb  = 128; % Nombre d'echantillons par symbole
D = 100; % Débit binaire en b/s
Tb = 1/D;   % Temps de bit
fs  = D*Nech_symb; % Frequence d'echantillonnage
Ts = 1/fs;  % Temps d'échantilonnage
EbNo = 1.6; % Rapport SNR en dB
Noise_Flag = 1;  

% Paramètres du filtre de mise en forme et génération de
% l'impulsion de Nyquist
alpha = 0.22;   % Rolloff
nyq_blocks = 8; % Duree de la réponse impulsionnelle

Ntot_ech = nyq_blocks*Nech_symb;
beta = alpha/(2*Tb);

% Variable temporelles
time_index=Ts*[-(Ntot_ech/2):(Ntot_ech/2 - 1)];
t=-nyq_blocks/2:1/Nech_symb:nyq_blocks/2-1/Nech_symb;  % Vecteur "temps"

% Evite les singularités
trouble_index = (abs(4*beta*time_index) == 1);
ok_index = ~trouble_index;

% Génération de l'impulsion
pulse = ones(size(time_index));
pulse(ok_index) = (cos(2*pi*beta*time_index(ok_index))./(1-(4*beta*time_index(ok_index)).^2));
g = sinc(time_index/Tb).*pulse;
Es=(sum(g.^2))/length(g);

% %%  PLOT de l'impulsion - Figure #1
plot(time_index, g);
grid

pause
%  Génération de la séquence d'information
M=2;      % Nombre d'états de la PAM
Lseq=1000;  % Longueur de la séquence d'information
 
if(M==2)
	M_ary_seq=2*floor(rand(1,Lseq)-0.5)+1; % Alphabet binaire
else
	M_ary_seq=2*(floor(M*rand(1,Lseq)))-M+1; % Alphabet M-aire
end

%  Generation des impulsions de Dirac
tt=0:1/Nech_symb:Lseq-1/Nech_symb; 
s = zeros(1,Lseq*Nech_symb);
for i = 0:Lseq-1
 s(1+(i*Nech_symb)) = M_ary_seq(i+1);
end

% %%  PLOT de la séquence d'information M-aire(10 symboles). Figure #2
figure
stem(tt(1:10*Nech_symb),s(1:10*Nech_symb));
pause

%  Génération de la séquence codée
S = [conv(g,s) 0];
S = S(1:length(s)); % On raccourci le transitoire final
if Noise_Flag
    N=2*sqrt(Es/(2*(10^(EbNo/10))))*randn(1,length(S)); % Bruit gaussien aditif
    S = S + N;
else
end
MM=max(max(S),abs(min(S)))+0.1; %0.1 de marge par rapport à l'amplitude max

% %%  PLOT de la séquence modulée.  Figure #3
Nsymb=20;  %nombre de symboles à représenter
temps=-nyq_blocks/2:1/Nech_symb:Lseq;
temps=temps(1:length(S));
plot(temps(1:Nsymb*Nech_symb),S(1:Nsymb*Nech_symb)) % les 20 premiers symboles
axis([-nyq_blocks/2 Nsymb-nyq_blocks -MM MM]);
set(gca,'XTick',-nyq_blocks/2:1:Nsymb-nyq_blocks);
grid
for i=1:2:M
 X=[temps(1) temps(Nsymb*Nech_symb)];
 Y=[i i];
 Ym=[(-1)*i (-1)*i];
 plot(X,Y,'r',X,Ym,'r')
end
pause

% Construction du DO
S = S(Ntot_ech+1:length(S)-Ntot_ech); %Troncature du transitoire
T = 0:1:length(S)-1;
Mult=length(S)/Nech_symb;
Eye_periods = 2;

SS = reshape(S,Eye_periods*Nech_symb,Mult/Eye_periods);
[lines,cols]=size(SS);
SSt=[SS(1,2:cols) SS(lines,cols)]; % 1ère ligne décalée et adjointe
SSS=[SS;SSt];
dt = Eye_periods/(2*Nech_symb);
tt = -1:dt:Eye_periods/2;

% %%  PLOT le DO
figure
plot(tt,SSS)
axis([-1 1 -MM MM])
grid
X=[0 0];
Y=[-MM MM];
plot(X,Y,'g');

