unipolar return to zero line code psd not equal to theoretical calculation

Hi, all

Now I am try to calculate unipolar RZ line code PSD by MATLAB, attached my code as following, for comparision, I input theoretical equation is

where

is variance of random input bits, is the means of random input bits.

but result is not equal, FFT result harmonic spectrum is not equal to the theoretical calculation, could you help to check my code? thanks.

clc; close all; clear all;

% waveform setting
ga=1;           % g(t) magnibude
gr=2;           % symbol rates,Mhz
gt=1/gr;        % symbol interval
gd=0.5;         %
gn=2^7;         % g(t) symbol number

gs=64;          % sampling points per symbol
gp=gn*gs;        % total sampling point

dt=gt/gs;        % symbo time interval
st=gp*dt;		% total sampling time

df=1/st;        % samping frequency interval

bw=gp*df;		% system bandwidth
t=linspace(-st/2,st/2,gp);	% time vector
f=linspace(-bw/2,bw/2,gp);        % freq vector

loop=1000;            % for smooth the curve
EP=zeros(1,gp);
for jj=1:loop
    a=round(rand(1,gn));    % generate gn symbols
    ma=mean(a,2);
    ta=std(a,0,2);          % variance
    [s t] = lineencoder('unirz',a,gt,ga,gs,0.5);
    O=fftshift(fft(s));
    Q=abs(O./gp);		      % time to freq, two side spectral 

    P=Q.^2;       %P
    EP=(EP*(jj-1)+P)/jj;  %????

end

% f1=figure('Name','pam_spectral_nrz_unipolar_f1','Visible','off');
% plot(t(1:10/gt*gs),s(1:10/gt*gs),'b');
% axis([1 10 -0.5 1.5]);
% xlabel('t (Sec)');ylabel('Volts(V)');title('RZ Unipolar Symbols');



f2=figure('Name','pam_spectral_nrz_unipolar_f2','Visible','on');
ps=30+10*log10(EP./df+eps);% plus eps to avoid device by zero, unit dBm so puls 30
plot(f,ps,'r');hold on;grid on;
xlabel('f (MHZ)');ylabel('Ps(f)');title('RZ Unipolar PSD by MATLAB FFT');
axis([-bw/8 bw/8 -60 60]);


f3=figure('Name','pam_spectral_nrz_unipolar_f3','Visible','on');

gf=ga*gt/2*sinc(f*gt/2).*exp(-1i*pi.*f*gt/2);
m=-bw/2:gr:bw/2-gr;
m=m/2;
gmf=zeros(1,bw/df);
% gmf0=sinc(m/2);
gmf(1:gr/df:end)=(ga*gt)./2*sinc(m/2);
p=ta^2/gt*abs(gf).^2 + (ma/gt)^2.*abs(gmf).^2;
% p=(ma/gt)^2.*abs(gmf0).^2;
% p=(ma/gt)^2.*abs(gmf).^2;


plot(f,10*log10(p)+30,'b');grid on;
axis([-bw/8 bw/8 -60 60]);




xlabel('f (MHZ)');ylabel('Ps(f)');title('RZ Unipolar PSD by Math Calculation');






function [x T] = LineEncoder(type,inbits,Tb,A,Fs,D)

if nargin<4, A = 1;end
if nargin<3, Tb = 1e-9;end
if nargin<2, inbits = [1 0 1 0];end
if nargin<1, type = 'uninrz';end
%---Implementation starts here
Rb = 1/Tb; %---Bit rate
% Fs = 4*Rb;
N = length(inbits);   %---Bit Length of input bits
tTb = linspace(0,Tb,Fs); %---interval of bit time period
x = [];
switch lower(type)
    case 'uninrz'
        for k = 1:N
            x = [x A*inbits(k)*ones(1,length(tTb))];
        end
    case 'unirz'
        for k = 1:N
            x = [x A*inbits(k)*ones(1,floor(length(tTb)*D)) 0*inbits(k)*ones(1,length(tTb)-floor(length(tTb)*D))];
        end
    case 'polrz'
        for k = 1:N
            c = ones(1,floor(length(tTb)*D));
            b = zeros(1,length(tTb)-floor(length(tTb)*D));
            p = [c b];
            x = [x ((-1)^(inbits(k)+1))*(A)*p];
        end
    case 'polnrz'
        for k = 1:N
            x = [x ((-1)^(inbits(k) + 1))*A*ones(1,length(tTb))];
        end
    case 'manchester'
        for k = 1:N
            c = ones(1,length(tTb)/2);
            b = -1*ones(1,length(tTb)/2);
            p = [c b];
            x = [x ((-1)^(inbits(k)+1))*A*p];
        end
    case 'ami'
end
T = linspace(0,N*Tb,length(x)); %---Time vector for n bits

end