a=[1/64 1/36 1/16 1/4];c=pi^2;
t=0:.025:1;
t0=clock;
for q=1:4
z=zeros(size(t)); s1=zeros(size(t));s2=zeros(size(t));
s3=zeros(size(t));s=zeros(size(t));s11=zeros(size(t));

for p=1:length(t)
for l=1:100
s1(p)=s1(p)+cos(2*pi*l*t(p))*exp(-2*(l*pi)^2*a(q));
end
end
for p=1:length(t),
for m=1:100,
for n=1:100, 
s2(p)=s2(p)+(1/(m*n))*(cos(2*pi*(m-n)*t(p))*exp(-2*(m-n)^2*c*a(q))-cos(2*pi*(m+n)*t(p))*exp(-2*(m+n)^2*c*a(q)));
end
end
end
for p=1:length(t),
for m1=1:100
for n1=1:100
s3(p)=s3(p)+(1/(m1*n1))*exp(-(m1^2+n1^2)*c*a(q))*sin(2*pi*n1*t(p))*sin(2*pi*m1*t(p));
end
end
end
for p=1:length(t),
s11(p)=(s1(p)*4+1)*a(q)+(1/c)*(.5*s2(p)-s3(p));

%end

a1(p)=10*log10(s11(p)*2/(2831.31^2))+87.18;
%m=mean(a1)
%sigma=std(a1)

y1=['s11',num2str(q),'(p)'];

y2=['=(s1(p)*4+1)*a(q)+(1/c)*(.5*s2(p)-s3(p));'];
y=[y1 y2];
eval(y)
v1=['a1',num2str(q),'(p)'];
v2=['=a1(p);'];
v=[v1 v2];
eval(v)
end
z=['  semilogy(t,s11',num2str(q),')'];
eval(z)

hold on
end
durata=etime(clock,t0)
xlabel('deriva, X/V');ylabel('Puterea la iesire, unitatea este V la patrat')
title('Variatia puterii zgomotului la iesirea cuantificatorului-legea miu- in functie de deriva normata')