si.m

很多matlab的源代码

function Y = si(x,a)
% SI    The sine integral.
%
%	Y = SI(X,A) Computes the integral of sin(At)/(At) over [0, X] 
%
%    	SI (with no input arguments) invokes the following example:
%
%	% Plot sinc(t)=sin(pi*t)/(pi*t) and the area under it over [0 5]
%    	  >>t = 0:0.05:5;
%    	  >>x1 = sinc(t);
%    	  >>x2 = si(t,pi);
%    	  >>plot(t,x1,t,x2)


% ADSP Toolbox: Version 2.0 
% For use with "Analog and Digital Signal Processing", 2nd Ed.
% Published by PWS Publishing Co.
%
% Ashok Ambardar, EE Dept. MTU, Houghton, MI 49931, USA
% http://www.ee.mtu/faculty/akambard.html
% e-mail: [e-mail]
% Copyright (c) 1998



if nargin==0,help si,disp('Strike a key to see results of the example')
pause,t=0:.05:5;x1=sinc(t);x2=si(t,pi);
 vx=matverch;
 if vx < 4, eval('clg');else,eval('clf');end
plot(t,x1,t,x2),axesn,grid
return,end

tol=2e-16;
if nargin<2,a=1;end
Y=0*x;s=prod(size(x));y=abs(x*a);z=y.*y;
i1=find(y<20);s1=prod(size(i1));
if s1>0
z1=z(i1);y1=y(i1);
Y(i1)=y1.*ghgf(1,2,[0.5 1 1 1],[1.5 1.5 1 1],-z1/4,tol);
end
if s>s1
i2=(1:s);i2(i1)=[];z1=z(i2);y1=y(i2);
Y1=ghgf(3,0,[1.5 1 1 1],[1 1 1 1],(-4)./z1,tol);
Y2=ghgf(3,0,[0.5 1 1 1],[1 1 1 1],(-4)./z1,tol);
Y3=0.5*pi-(sin(rem(y1,2*pi)).*Y1./y1 + cos(rem(y1,2*pi)).*Y2)./y1;Y(i2)=Y3;
end
Y=Y.*(y~=0).*sign(x)/abs(a);