DSP for dispersion compensation by FFT overlap-add method

I have tried to achieve a simulation of 100G DP-QPSK system. This simulation uses FFT overlap-add method to compensate dispersion.The 100G DP-QPSK system was simulated by Optisystem sofware. DSP for compensate dispersion used matlap code.

But it was not successful!

This is my matlap code.

(main.m)

close all; format short;

OutputPort1 = InputPort1; OutputPort2 = InputPort2; OutputPort3 = InputPort3; OutputPort4 = InputPort4;

%%%%%%%%%%%%%%%%%%%%%%% Parameters Definition %%%%%%%%%%%%%%%%%%%%%%%%%%%%% ChannelWavelength = ChannelWavelength * 1e-9; DispersionRefWavelength = DispersionRefWavelength * 1e-9; ResidualDispersion = ResidualDispersion * 1e-12 / 1e-6; ResidualDispersionSlope = ResidualDispersionSlope * 1e-12 / 1e-15; c0 = 2.99792458e8; fft_length=128;

%%%%%%%%%%%%%%%%%%%%%%%%% Down Sampling %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% SamplesPerSymbolIn = length(InputPort1.Sampled.Signal) / NumberOfSymbols; % calculates number of samples per symbol

InX_Signal = ones(1, NumberOfSymbols * SamplesPerSymbol); InX_Noise = ones(1, NumberOfSymbols * SamplesPerSymbol); InY_Signal = ones(1, NumberOfSymbols * SamplesPerSymbol); InY_Noise = ones(1, NumberOfSymbols * SamplesPerSymbol);

for m = 1 : NumberOfSymbols * SamplesPerSymbol

    InX_Signal(1,m) =   InputPort1.Sampled.Signal((2*m-1) * SamplesPerSymbolIn / (2*SamplesPerSymbol)) + 1i * InputPort2.Sampled.Signal((2*m-1) * SamplesPerSymbolIn / (2*SamplesPerSymbol));
    InX_Noise(1,m)  =   InputPort1.Noise.Signal((2*m-1) * SamplesPerSymbolIn / (2*SamplesPerSymbol)) + 1i * InputPort2.Noise.Signal((2*m-1) * SamplesPerSymbolIn / (2*SamplesPerSymbol));

    InX             =   InX_Signal + InX_Noise;

    InY_Signal(1,m) =   InputPort3.Sampled.Signal((2*m-1)*SamplesPerSymbolIn/(2*SamplesPerSymbol))+1i*InputPort4.Sampled.Signal((2*m-1)*SamplesPerSymbolIn/(2*SamplesPerSymbol));
    InY_Noise(1,m)  =   InputPort3.Noise.Signal((2*m-1)*SamplesPerSymbolIn/(2*SamplesPerSymbol))+1i*InputPort4.Noise.Signal((2*m-1)*SamplesPerSymbolIn/(2*SamplesPerSymbol));

    InY             =   InY_Signal + InY_Noise;
end;

clear InX_Signal; clear InX_Noise; clear InY_Signal; clear InY_Noise;

%In_X: signal and noise for polarization X %In_Y: signal and noise for polarization Y

SampledVector = (2*(1:NumberOfSymbols*SamplesPerSymbol)-1)*SamplesPerSymbolIn/(2*SamplesPerSymbol); TimeSampledVector = InputPort1.Sampled.Time(SampledVector);

%%%%%%%%%%%%%% Normalizing %%%%%%%%%%%%%%%%%%%%%% OutX = InX/sqrt(mean(abs(InX).^2)); OutY = InY/sqrt(mean(abs(InY).^2));

StoreInX = OutX; StoreInY = OutY;

clear SampledVector; %%%%%%%%%%%%%%%%%%%%%% Dispersion Compensation in Frequency domain %%%%%%%%%%%%%%%

if EnableD == 1 data_IQx=OutX; data_IQy=OutY; %sample_period=1e12/(SymbolRate*NewSamplesPerSymbol); %recovering sampling period in 'ps' sample_period=1e12/(SymbolRate*SamplesPerSymbol); Dispersion=ResidualDispersion*PropagationLength; %CD_fil_fde generates the filter and perform the filtering in the %frequency domain [data_IQx data_IQy]=CD_fil_fde (data_IQx, data_IQy, fft_length, Dispersion, DispersionRefWavelength*1e9, sample_period);

%%Nonlinear compensator NLC = 0; if NLC == 1

        alpha= 25; 
        beta = 25;
        %they need to be optimized
        data_IQx=data_IQx.*exp(-1i*(alpha*(abs(data_IQx)).^2+beta*(abs(data_IQy)).^2));
        data_IQy=data_IQy.*exp(-1i*(alpha*(abs(data_IQy)).^2+beta*(abs(data_IQx)).^2));
    end

end OutX=data_IQx; OutY=data_IQy; AfterDX = OutX; AfterDY = OutY; %%%%%%%%%%%%%%%%%%%%%%%%%%%%% Interpolation and Output to OptiSystem %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% OutputPort1.Sampled.Signal = interp1(TimeSampledVector,real(OutX),InputPort1.Sampled.Time,InterpolationMethod); OutputPort2.Sampled.Signal = interp1(TimeSampledVector,imag(OutX),InputPort1.Sampled.Time,InterpolationMethod); OutputPort3.Sampled.Signal = interp1(TimeSampledVector,real(OutY),InputPort1.Sampled.Time,InterpolationMethod); OutputPort4.Sampled.Signal = interp1(TimeSampledVector,imag(OutY),InputPort1.Sampled.Time,InterpolationMethod);

%% OVERLAP-ADD METHOD %CD_fil_fde CD compensation in frequency domain . %________Algorithmic details______________ %This function construct the frequency filter and call the Overlap-add %method main function of filtering using FFT. %datax , datay are the input data stream . %fft_length is the length of the FFT, L is the length of the data window %that we consider and G is the total zero pad length that is applied before %and after the data window . function [ data_IQx data_IQy ]=CD_fil_fde ( datax , datay , fft_length , dispersion, wavelength , sample_period ) m = size ( datax , 1) ; if m == 1 datax = datax ( : ) ; datay = datay ( : ) ; end m = size (datax , 1) ; if m == 1 datay = datay ( : ) ; end

ndata = size ( datax , 1 ) ;

if fft_length > 2^20 , err ( generatemsgid ( ' filterTooLong ' ) , ' Filters of length greater than 2^20 are not supported . Use dfilt .fftfir instead . ' ) ; end if fft_length > ndata err (generatemsgid ( 'filterTooShort ' ) , ' Filter length must be shorter than data length ' ) ; else %% L and G calculation

G = fft_length / 4 ; L =fft_length / 2 ; %% construction of the filter %here the filter is constructed using the frequency domain formula directly %from the theory

C = 3e5 ; %speed of light in km/s q=(-fft_length / 2 : fft_length /2 - 1)' ; omega=2*pi*(q) /( sample_period * fft_length ) ; filt=exp(-1i ( dispersion *( wavelength ) ^2/(4 pi *C) ) *omega .^ 2 ) ;

%fil is a structure where are saved all the parameters needed for perform the o v e r la p method fil = struct ( 'filter' , filt , 'nfft' , fft_length , 'window' , L, 'ndata' , ndata ,'overlap' , G) ;

%main program data_IQx = overlap_add ( fil , datax ) ; data_IQy = overlap_add ( fil , datay ) ; end end

%% Overlap_add method of FIR filtering using FFT. % -------Algorithmic details---------------- % b is a structure where are stored all the parameter needed to perform the overlap_add method . % b = struct ( ' filter ' , . . . , ' nfft ' , . . . , 'window ' , . . . , ' ndata ' , . . . , 'over lap ' , . . . ) % filter is the filter defined in frequency domain , nfft is the FFT % length , window is the window data length , ndata is the length of all % the received data , overlap is the zero padding length defined as nfft/4 % The overlap/add algorithm convolves b.filter with blocks of data , and adds % the overlapping output blocks . I t uses the FFT to compute the % convolution . function y = overlap_add (b , data ) %% l o a d ing parameters

nfft=b.nfft ; ndata=b.ndata ; L=b.window ; G=b.overlap ; B = fftshift (b.filter) ; if length(B)==1, B = B(:) ; % make sure that B is a column end if size (B,2)==1 B = B(:,ones ( 1 , size ( data , 2 ) ) ) ; % replicate the column B end if size ( data , 2 )==1 data = data ( : , ones ( 1 , size (b , 2 ) ) ) ; % replicate the column data end y = zeros ( size ( data ) ) ; istart = 1 ; while istart <= ndata iend = min ( istart+L-1,ndata ) ; if ( iend - istart ) <= 0 break else X = fft([(zeros ( 1 ,G))'; data( istart : iend , : )] , nfft ) ; end Y = ifft (X.*B) ; yend = min( ndata , istart+nfft -1) ; y ( istart : yend , : ) = y ( istart : yend , : ) + Y( 1 : ( yend-istart +1) , : ) ; istart = istart + L; end end