tic
j=2048*52; % this is the frequency about which we are summing
n=25; % this is the frequency width to the left and right of r
N=4096*52;
w=[1/2;ones(25,1);1/2];     
T=[];
        for alpha=-1:1
             for beta=-1:1
                 for alphap=-1:1
                    for betap=-1:1
                        for k=(j-13):(j+13)
                            for kp=(j-13):(j+13)  
PhiaS=1/4*(-1/2)^(abs(alpha)+abs(alphap))*sin(pi/2*(k-kp+2*alpha-2*alphap))*exp(i*pi*(alpha-alphap)*(1-1/N));
Phib=1/4*(-1/2)^(abs(beta)+abs(betap))*sin(pi/2*(k-kp+2*beta-2*betap))*exp(-i*pi*(beta-betap)*(1-1/N));
Den1=sin(pi/(2*N)*(k-kp+2*alpha-2*alphap));
Den2=sin(pi/(2*N)*(k-kp+2*beta-2*betap));
if k-kp+2*alpha-2*alphap==0
R1=1/4*(-1/2)^(abs(alpha)+abs(alphap))*N*exp(i*pi*(alpha-alphap)*(1-1/N));
else R1=PhiaS/Den1;
end
if k-kp+2*beta-2*betap==0
R2=1/4*(-1/2)^(abs(beta)+abs(betap))*N*exp(-i*pi*(beta-betap)*(1-1/N));
else R2=Phib/Den2;
end
T=[T w(k-(j-13-1))*w(kp-(j-13-1))*R1*R2];
                             end
                         end
                     end  
                 end   
             end
         end
ZZ=sum(T)/N^2;
toc