rat.m

可用于对信号的频谱进行分析,希望站长能整理好,以便大家互相学习

function [N,D] = rat(X,tol)
%RAT    Rational approximation.
%   [N,D] = RAT(X,tol) returns two integer matrices so that N./D
%   is close to X in the sense that abs(N./D - X) <= tol*abs(X).
%   The rational approximations are generated by truncating continued
%   fraction expansions.   tol = 1.e-6*norm(X(:),1) is the default.
%
%   S = RAT(X) or RAT(X,tol) returns the continued fraction 
%   representation as a string.
%
%   The same algorithm, with the default tol, is used internally
%   by MATLAB for FORMAT RAT.
%
%   Class support for input X:
%      float: double, single
%
%   See also FORMAT, RATS.

%   Copyright 1984-2004 The MathWorks, Inc. 
%   $Revision: 5.18.4.4 $  $Date: 2004/07/05 17:02:10 $

% Approximate x by
%
%                              1
%         d1 + ----------------------------------
%                                 1
%              d2 + -----------------------------
%                                   1
%                   d3 + ------------------------
%                                      1
%                        d4 + -------------------
%                                        1
%                             d5 + --------------
%                                           1
%                                  d6 + ---------
%                                             1
%                                       d7 + ----
%                                             d8
%
% The d's are obtained by repeatedly picking off the integer part and 
% then taking the reciprocal of the fractional part.  The accuracy of
% the approximation increases exponentially with the number of terms
% and is worst when x = sqrt(2).  For x = sqrt(2), the error with k
% terms is about 2.68*(.173)^k, which is
%
%         1    4.6364e-01
%         2    8.0210e-02
%         3    1.3876e-02
%         4    2.4006e-03
%         5    4.1530e-04
%         6    7.1847e-05
%         7    1.2430e-05
%         8    2.1503e-06
%         9    3.7201e-07
%        10    6.4357e-08

if nargin < 2
    tol = 1.e-6*norm(X(isfinite(X)),1);
end
if ~isreal(X)
    if norm(imag(X),1) <= tol*norm(real(X),1)
        X = real(X);
    elseif nargout > 1
        [NR,DR] = rat(real(X));
        [NI,DI] = rat(imag(X));
        D = lcm(DR,DI);
        N = D./DR.*NR + D./DI.*NI*i;
        return
    else
        N = strvcat(rat(real(X)),' +i* ...',rat(imag(X)));
        return
    end
end
if nargout > 1
    N = zeros(size(X), class(X));
    D = zeros(size(X), class(X));
else
    N = zeros(0,0,class(X)); 
end
for j = 1:numel(X)
    x = X(j);
    if ~isfinite(x), % Special case for inf, -inf, NaN
        if nargout <=1
            s = int2str(x);
            k = length(s)-size(N,2);
            N = char([N ' '*ones(j-1,k); s ' '*ones(1,-k)]);
        else
            if ~isnan(x), N(j) = sign(x); else N(j) = 0; end
            D(j) = 0;
        end
    else
        k = 0;
        C = [1 0; 0 1];    % [n(k) n(k-1); d(k) d(k-1)];
        while 1
            k = k + 1;
            neg = x < 0;
            d = round(x);
            if ~isinf(x),
                x = x - d;
                C = [C*[d;1] C(:,1)];
            else % Special case for +/- inf
                C = [[x;0] C(:,1)];
            end
            if nargout <= 1
                d = int2str(abs(d));
                if neg, d = ['-' d]; end
                if k == 1
                    s = d;
                elseif k == 2
                    s = [s ' + 1/(' d ')'];
                else
                    p = min(find(s==')'));
                    s = [s(1:p-1) ' + 1/(' d ')' s(p:p+k-3)];
                end
            end
            if (x==0) | (abs(C(1,1)/C(2,1) - X(j)) <= max(tol,eps(X(j))))
               break
            end
            x = 1/x;
        end
        if nargout > 1
            N(j) = C(1,1)/sign(C(2,1));
            D(j) = abs(C(2,1));
        else
            k = length(s)-size(N,2);
            N = char([N ' '*ones(j-1,k); s ' '*ones(1,-k)]);
        end
    end
end