% Example of numerical simulation. System: fiber with only GVD
%
%                   __
%                  /  \                      
%                 |    |       
%                  \  /       
%      x(t)  ----------------  y(t)
%                  fiber       
%                                                  
%   Author: Paolo Serena, Mar-2006
%   University of Parma, Italy

clear all
%close all   % close all figures
clc         % clear screen

%%% Fiber data. They are normalized to a reference time T equal to the 1/e
% width. If p(t) is the power of a Gaussian pulse, t being the non-normalized  
% time, it is p(T/2)=1/exp(1). In normalized units it is p(0.5)=1/exp(1).

zoverLd = 0.5;        % z/Ld
signb2 = 1;         % sign(beta2)

%%% Numerical data
Nt = 32;    % number of points x bit (points or samples between t=0 and t=T)
Nbit = 16;  % number of bits. Note: Total Number of discrete points: Nt*Nbit
m = 1;      % super Gaussian order.
%%%%%%%%%%%%%%%%%%%%%%%%%% program %%%%%%%%%%%%%%%%%%%%%%%%%%

fn = fftshift(-Nt/2:1/Nbit:Nt/2-1/Nbit);    % frequencies normalized to the 
                                            % the reference frequency R=1/T. 
                                            
% Note: 1) Lowest frequency (resolution) = 1/Nbit        
%       2) Largest frequency (half of Nyquist frequency) = Nt/2

time = -Nbit/2:1/Nt:Nbit/2-1/Nt;    % time normalized to 1/R and 
                                    % centered around 0 (useful convention)
                            
x = exp(-(2*time).^(2*m)/2);    % signal: Gaussian pulse (centered around 0)
% Note: different convention than Agrawal ("Nonlinear Fiber Optics", AP).
% Thanks to the first factor 2, here we can call "bit time" the reference
% time T. For this reason a Gaussian pulse doubles its width at zoverLd=1/2.

xf = fft(x);            % frequency domain
wn = 2*pi*fn;                                            
fiber = exp(-i*0.5*signb2*zoverLd*wn.^2);    % fiber transfer function
y = ifft(xf .* fiber);               % output signal in the time  domain

figure(1)
plot(fftshift(fn),20*log10(abs(fftshift(xf))))  % spectrum
xlabel('frequency   [a.u.]')
ylabel('|fft(x)|    [dB]')
grid on
hold on

figure(2)
subplot(2,1,1)
plot(time,abs(x).^2,'b',time,abs(y).^2,'r')
grid on
hold on
ylabel('Power  [a.u.]')
subplot(2,1,2)
plot(time,unwrap(angle(y)),'r') % unwrap to 2*pi phase jumps
grid on
hold on
ylabel('Phase  [a.u.]')
xlabel('time   [a.u.]')