invsi2.m

很多matlab的源代码

function [b,k] = invsi2(y,a,tol,kmax)
% INVSI2 Inverse sine square integral.
%
%	[C, K] = INVSI2(X,A,TOL,ITER) Inverse of the sine square integral SI2
%	X = SI2(y,A) where SI2 is integral of [sin(At)/(At)] over [0, y
%	A is the argument of At in sin(At)/(At) [Default: A = 1] 
%	TOL = desired tolerance [Default: TOL = 1e-7], 
%       ITER = maximum number of iterations [Default: ITER = 100]
%	C returns the inverse (which should equal X to within specifed TOL)
%       K returns the actual number of iterations performed.
%
%       INVSI2 (with no input arguments) invokes the following example:
%
%       % Find t at which 50%, 90%, and 95% of the maximum energy is contained
%	% in x(t) = sinc(t)=sin(pi*t)/(pi*t).
%	% Since E = area under x^2(t) over [-inf,  inf] = 1 (we know this), 
%	% E for [0 inf] = 0.5 and 50%, 90%, 95% of E equal 
%	% we use these to compute 
%	  >>e = [0.25 0.45 0.475]
%         >>t = invsi2([0.25 0.45 0.475],pi)
%       % As a check, e1=si2(t,pi) should equal e (to within TOL=1e-7) 
%          >>e1 = si2(t,pi)
%	   >>err = e-e1


% 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 invsi2,disp('Strike a key to see results of the example')
pause,e=[0.25 0.45 0.475],t=invsi2(e,pi),e1=si2(t,pi),err=e-e1,return,end

if nargin<4,kmax=100;end,
if nargin<3,tol=1e-7;end,
if nargin<2,a=1;end
fl1=0;fl2=0;x1=abs(a*y);x=x1.*(x1<pi/2);
[m,n]=size(x);bp=ones(m,n);b=0*bp;
k=0;
while any(any(abs(b-bp) > tol)) & k<=kmax
k=k+1;bp=b.*(b>=0)+.5*bp.*(b<0);
bp1=rem(bp,2*pi);f=(si2(bp)-x)/2;
fd=(sin(bp1+eps)./(bp+eps)).^2;b=(bp-f./fd);
end
if k>kmax,disp(['No convergence in ' num2str(kmax) ' iterations']),end
b=b.*sign(y)/abs(a);
if any(any(x1==pi/2)),i=find(x1==pi/2);b(i)=inf*sign(y(i));end
if any(any(x1>pi/2)),i=find(x1>pi/2);[p,l]=size(i);b(i)=NaN*ones(p,l);end