function [T] = oko(EbNo,B,PMD,L,typ);

%################# © Ing. Martin Kyselak, Bc. Petr Dorociak ###############
%# Skript pro zobrazení vlivu Polarizacni vidove disperze na pokles SNR pro
%ruzne druhy modulaci 
%# Funkci zadavejte ve tvaru [T] = oko(EbNo,B,PMD,L,typ), kde
%# EbNo = pomer bitove energie od sumove spectralni hustoty, napr.: [10:15]
%# B = prenosova rychlost [bit/s]
%# PMD = koeficient PMD [s/sqrt(km)]
%# L = delka vlakna [km]
%# typ = znaci typ modulace: 
%# 0 = NRZ-OOK (ASK) modulace
%# 1 = RZ-OOK (ASK) modulace
%# 2 = NRZ-DBPSK modulace
%# 3 = RZ-DBPSK modulace
%# 4 = NRZ-DQPSK modulace 
%# 5 = RZ-DQPSK modulace 
%################# © Ing. Martin Kyselak, Bc. Petr Dorociak ###############
 
help(mfilename);
close all hidden;
clc;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%    Maximalni delka vlakna pri zadane PMD a bitove rychlosti        %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
PMD_ps = PMD*1e12;
B_gb = B*1e-9;

%vypocet maximalni delky 
Lmax = 0.01./((B.*PMD).^2);

%vypocet DGD
DGD = PMD*sqrt(L);

%prevod na ps
DGD_ps = DGD*1e12;   

%generace osy x
x1 =0:PMD/100:PMD;

% vypocet logaritmicke krivky
Lmax_log = 0.01./((B.*x1).^2);

%vykresleni logaritmicke krivky
semilogy(x1, Lmax_log, '-','LineWidth',2);

grid on;
hold on;

if L <= Lmax 

%vykresleni bodu maximalni delky
semilogy(PMD, Lmax, '*r','LineWidth',2);
title( 'Zavislost maximalni delky na PMD a prenosove rychlosti');
xlabel( '\rightarrow  {\it PMD} [s*sqrt(km)]');
ylabel( '\rightarrow {\it Lmax} [km]'); 
legend('Logaritmicka krivka',['Maximalni delka vlakna Lmax = ' num2str(Lmax) ' km, pri B = ' num2str(B_gb) ' Gb/s']);  

% vypis na displej
disp(['Maximalni delka vlakna pro zadanou PMD = ' num2str(PMD_ps) ' ps/sqrt(L) a bitovou rychlost B = ' num2str(B_gb) ' Gb/s je Lmax = ' num2str(Lmax) ' km.']);
disp('  ');  

else
    
%vypocet maximalni PMD
PMDmax = sqrt(0.01./((B.^2.*L)));

semilogy(PMD, Lmax, '*r','LineWidth',2);
semilogy(PMDmax, L, '*g','LineWidth',2);
title( 'Zavislost maximalni delky na PMD a prenosove rychlosti');
xlabel( '\rightarrow  {\it PMD} [s*sqrt(km)]');
ylabel( '\rightarrow {\it Lmax} [km]'); 
legend('Logaritmicka krivka',['Maximalni delka vlakna Lmax = ' num2str(Lmax) ' km, pri B = ' num2str(B_gb) ' Gb/s'], ['Minimalni hodnota PMD pri zadane L = ' num2str(L) ' km, pri B = ' num2str(B_gb) ' Gb/s'] );  

% vypis na displej
disp(['Maximalni delka vlakna pro zadanou PMD = ' num2str(PMD_ps) ' ps/sqrt(L) a bitovou rychlost B = ' num2str(B_gb) ' Gb/s je Lmax = ' num2str(Lmax) ' km.']);
disp('  ');  
disp(['Vami zadana delka L = ' num2str(L) ' km je vetsi nez maximalni mozna delka Lmax = ' num2str(Lmax) ' km vlakna pro zadanou PMD = ' num2str(PMD_ps) ' ps/sqrt(L) a bitovou rychlost B = ' num2str(B_gb) ' Gb/s']);
disp('  ');  

end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%                     NRZ-OOK (ASK) modulace                         %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% typ:OOK (ASK) modulace  
if typ == 0  

%pocet stavu
M=2;

%vypocet BER
ber = berawgn(EbNo,'pam',M);

% konvertovani na SNR - Poznamka, protoze No = 2noiseVariance^2, pridame k SNR pro PAM unipolarni 6 dB 
SNR = EbNo + 6 + 10*log10(log2(M));

% vykresleni grafu
figure(2);
semilogy(SNR,ber,'b','LineWidth',2);
title( 'SNR v zavislosti na BER u NRZ-OOK (ASK) modulace');
xlabel( '\rightarrow {\it SNR} [dB]');
ylabel( '\rightarrow {\it BER} [-]');

grid on; 
drawnow;
hold on;   

% vypocitani SNR penalty
g=0.5; %power splitting ratio
A=22; %Raissed cossine filter

T=1/B;%bit period = Tsp

normDGD=DGD/T;%normalized DGD

%vypocet Eye penalty
P_ook_nrz = 2*A*(normDGD)^2*(g*(1-g));

%SNR penalty
SNRP=SNR-P_ook_nrz;

% vykresleni grafu
semilogy(SNRP,ber,'r','LineWidth',2);

legend('SNR pri DGD = 0 ps' , ['SNR pri DGD = ' num2str(DGD_ps) ' ps']); 

grid on; 
drawnow;
hold on;

%vypis na displej
disp(['Power penalty pro OOK-NRZ = ' num2str(P_ook_nrz) ' dB pri DGD = ' num2str(DGD_ps) 'ps, respektive pri normDGD = ' num2str(normDGD) ' ']);


%otevreni modelu v simulinku
open('D:\PMD\ook_nrz.mdl');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%                      RZ-OOK (ASK) modulace                         %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% typ:OOK (ASK) modulace  
elseif typ == 1  

%pocet stavu
M=2;

%vypocet BER
ber = berawgn(EbNo,'pam',M);

% konvertovani na SNR - Poznamka, protoze No = 2noiseVariance^2, pridame k  SNR pro PAM unipolarni 6 dB 
SNR = EbNo + 6 + 10*log10(log2(M));

% vykresleni grafu
figure(2);
semilogy(SNR,ber,'b','LineWidth',2);
title( 'SNR v zavislosti na BER u RZ-OOK (ASK) modulace');
xlabel( '\rightarrow {\it SNR} [dB]');
ylabel( '\rightarrow {\it BER} [-]');

grid on; 
drawnow;
hold on;   

% vypocitani SNR penalty
g=0.5; %power splitting ratio
A=22; %Raissed cossine filter

T=1/B;%bit period = Tsp
normDGD=DGD/T;%normalized DG

%vypocet Eye penalty
P_ook_rz = (2*A*(normDGD)^2*(g*(1-g)))/2;

%SNR penalty
SNRP=SNR-P_ook_rz;


% vykresleni grafu

semilogy(SNRP,ber,'r','LineWidth',2);
legend('SNR pri DGD = 0 ps' , ['SNR pri DGD = ' num2str(DGD_ps) ' ps' ]); 

%vypis na displej
disp(['Power penalty pro OOK-RZ = ' num2str(P_ook_rz) ' dB pri DGD = ' num2str(DGD_ps) 'ps, respektive pri normDGD = ' num2str(normDGD) ' ']);

grid on; 
drawnow;
hold on;

%otevreni modelu v simulinku
open('D:\PMD\ook_rz.mdl');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%                     NRZ-DNPSK modulace                             %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% typ:DBPSK modulace    
elseif typ == 2  

%pocet stavu
M= 2; 

%vypocet BER
ber = berawgn(EbNo,'dpsk',M); 

% konvertovani na SNR
SNR = EbNo + 10.*log10(log2(M));

% vykresleni grafu
figure(2);
semilogy(SNR,ber,'b','LineWidth',2);
title( 'SNR v zavislosti na BER u NRZ-DBPSK modulace');
xlabel( '\rightarrow {\it SNR} [dB]');
ylabel( '\rightarrow {\it BER} [-]');
hold on;

% vypocitani SNR penalty
g=0.5; %power splitting ratio
A=25; %Raissed cossine filter

T=1./B;%bit period = Tsp
normDGD=DGD/T;%normalized DGD

%vypocet Eye penalty
P_dbpsk_nrz = 2.*A.*(normDGD).^2.*(g.*(1-g));

%SNR penalty
SNRP=SNR-P_dbpsk_nrz;


% vykresleni grafu

semilogy(SNRP,ber,'r','LineWidth',2);
legend('SNR pri DGD = 0 ps' , ['SNR pri DGD = ' num2str(DGD_ps) ' ps' ]); 

%vypis na displej
disp(['Power penalty pro DBPSK-NRZ = ' num2str(P_dbpsk_nrz) ' dB pri DGD = ' num2str(DGD_ps) 'ps, respektive pri normDGD = ' num2str(normDGD) ' ']);

grid on; 
drawnow;
hold on;

%otevreni modelu v simulinku
open('D:\PMD\dbpsk_nrz.mdl');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%                      RZ-DBPSK modulace                             %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% typ:DBPSK modulace    
elseif typ == 3  

%pocet stavu
M= 2; 

%vypocet BER
ber = berawgn(EbNo,'dpsk',M); 

% konvertovani na SNR
SNR = EbNo + 10.*log10(log2(M));

% vykresleni grafu
figure(2);
semilogy(SNR,ber,'b','LineWidth',2);
title( 'SNR v zavislosti na BER u RZ-DBPSK modulace');
xlabel( '\rightarrow {\it SNR} [dB]');
ylabel( '\rightarrow {\it BER} [-]');
hold on;

% vypocitani SNR penalty
g=0.5; %power splitting ratio
A=25; %Raissed cossine filter

T=1./B;%bit period = Tsp
normDGD=DGD/T;%normalized DGD

%vypocet Eye penalty
P_dbpsk_rz = (2.*A.*(normDGD).^2.*(g.*(1-g)))/2;

%SNR penalty
SNRP=SNR-P_dbpsk_rz;


% vykresleni grafu

semilogy(SNRP,ber,'r','LineWidth',2);
legend('SNR pri DGD = 0 ps' , ['SNR pri DGD = ' num2str(DGD_ps) ' ps' ]); 

%vypis na displej
disp(['Power penalty pro DBPSK-RZ = ' num2str(P_dbpsk_rz) ' dB pri DGD = ' num2str(DGD_ps) 'ps, respektive pri normDGD = '  num2str(normDGD) ' ']);

grid on; 
drawnow;
hold on;

%otevreni modelu v simulinku
open('D:\PMD\dbpsk_rz.mdl');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%                     NRZ-DQPSK modulace                             %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% typ:DQPSK modulace
elseif typ == 4 

%pocet stavu
M=4;

%vypocet BER
ber = berawgn(EbNo,'dpsk',M);

% konvertovani na SNR
SNR = EbNo + 10*log10(log2(M));

% vykresleni grafu
figure(2);
semilogy(SNR,ber,'b','LineWidth',2);
title( 'SNR v zavislosti na BER u NRZ-DQPSK modulace');
xlabel( '\rightarrow {\it SNR} [dB]');
ylabel( '\rightarrow {\it BER} [-]');

grid on; 
drawnow;
hold on;

% vypocitani SNR penalty
g=0.5; %power splitting ratio
A=25; %Raissed cossine filter

T=1/B;%bit period = Tsp
normDGD=DGD/T;%normalized DGD

%vypocet Eye penalty
P_dqpsk_nrz = (2*A*(normDGD)^2*(g*(1-g)))/4;

%SNR penalty
SNRP=SNR-P_dqpsk_nrz;


% vykresleni grafu

semilogy(SNRP,ber,'r','LineWidth',2);
legend('SNR pri DGD = 0 ps' , ['SNR pri DGD = ' num2str(DGD_ps) ' ps' ]); 

%vypis na displej
disp(['Power penalty pro DQPSK-NRZ = ' num2str(P_dqpsk_nrz) ' dB pri DGD = ' num2str(DGD_ps) 'ps, respektive pri normDGD = ' num2str(normDGD) ' ']);

grid on; 
drawnow;
hold on;

%otevreni modelu v simulinku
open('D:\PMD\dqpsk_nrz.mdl');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%                      RZ-DQPSK modulace                             %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% typ:DQPSK modulace
elseif typ == 5 

%pocet stavu
M=4;

%vypocet BER
ber = berawgn(EbNo,'dpsk',M);

% konvertovani na SNR
SNR = EbNo + 10*log10(log2(M));

% vykresleni grafu
figure(2);
semilogy(SNR,ber,'b','LineWidth',2);
title( 'SNR v zavislosti na BER u RZ-DQPSK modulace');
xlabel( '\rightarrow {\it SNR} [dB]');
ylabel( '\rightarrow {\it BER} [-]');

grid on; 
drawnow;
hold on;

% vypocitani SNR penalty
g=0.5; %power splitting ratio
A=25; %Raissed cossine filter

T=1/B;%bit period = Tsp
normDGD=DGD/T;%normalized DGD

%vypocet Eye penalty
P_dqpsk_rz = (2*A*(normDGD)^2*(g*(1-g)))/8;

%SNR penalty
SNRP=SNR-P_dqpsk_rz;


% vykresleni grafu

semilogy(SNRP,ber,'r','LineWidth',2);
legend('SNR pri DGD = 0 ps' , ['SNR pri DGD = ' num2str(DGD_ps) ' ps' ]); 

%vypis na displej
disp(['Power penalty pro DQPSK-RZ = ' num2str(P_dqpsk_rz) ' dB pri DGD = ' num2str(DGD_ps) 'ps, respektive pri normDGD = ' num2str(normDGD) ' ']);

grid on; 
drawnow;
hold on;

open('D:\PMD\dqpsk_rz.mdl');

end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%      Zavislost Power penalty na DGD/T pro ruzne druhy modulaci     %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
normDGD_x = 1;
g=0.5; %power splitting ratio

%generace osy x
x =0:normDGD_x/100:normDGD_x;

% OOK konstanty filtru
A_ook=22; %Raissed cossine filter

% PSK konstanta filtu
A_psk=25; %Raissed cossine filter

%vypocty Pe
Pe_ook_nrz = 2.*A_ook.*(x).^2.*(g.*(1-g));
Pe_ook_rz = (2.*A_ook.*(x).^2.*(g.*(1-g)))/2;
Pe_dbpsk_nrz = 2.*A_psk.*(x).^2.*(g.*(1-g));
Pe_dbpsk_rz = (2.*A_psk.*(x).^2.*(g.*(1-g)))/2;
Pe_dqpsk_nrz = (2.*A_psk.*(x).^2.*(g.*(1-g)))/4;
Pe_dqpsk_rz = (2.*A_psk.*(x).^2.*(g.*(1-g)))/8;

figure(3);
plot(x,Pe_ook_nrz,'-r','LineWidth',2);
hold on;
plot(x,Pe_ook_rz,'--r','LineWidth',2);
hold on;
plot(x,Pe_dbpsk_nrz,'-b','LineWidth',2);
hold on;
plot(x,Pe_dbpsk_rz,'--b','LineWidth',2);
hold on;
plot(x,Pe_dqpsk_nrz,'-m','LineWidth',2);
hold on;
plot(x,Pe_dqpsk_rz,'--m','LineWidth',2);
hold on;
grid on;
title( 'Zavislost Power penalty na DGD/Tsp pro ruzne druhy modulaci');
xlabel( '\rightarrow {\it DGD/Tsp} [-]');
ylabel( '\rightarrow {\it Power penalty} [dB]');
legend('OOK-NRZ','OOK-RZ','DBPSK-NRZ','DBPSK-RZ','DQPSK-NRZ','DQPSK-RZ'); 

