www.gusucode.com > signal 工具箱matlab源码程序 > signal/intfilt.m

    function [h,a]=intfilt(R,L,freqmult)
%INTFILT Interpolation FIR Filter Design.
%   B = INTFILT(L,P,ALPHA) designs a linear phase FIR filter which, when 
%   used on a sequence interspersed with L-1 consecutive zeros every L
%   samples, performs bandlimited interpolation using the nearest 2*P
%   non-zero samples and assuming an original bandlimitedness of ALPHA
%   times the Nyquist frequency.  B is length 2*L*P-1.
%
%   The bandlimitedness factor expresses how much of the Nyquist interval
%   is occupied by the signal to be interpolated. This factor should be
%   greater than zero and less or equal to one. A bandlimitedness factor
%   that is less than one allows for a larger transition band and for
%   "don't-care" regions in the stopband of the filter where the signal has
%   no significant spectral content. The net result of specifying ALPHA
%   less than one is better stopband attenuation for a given L and P. See
%   example below.
%
%   B = INTFILT(L,N,'Lagrange') designs an FIR filter which performs Nth 
%   order Lagrange polynomial interpolation on a sequence that is inter-
%   spersed with L-1 consecutive zeros every L samples.  B is length 
%   (N+1)*L-1 for N odd and (N+1)*L for N even.  If both N and L are even, 
%   the filter designed is NOT linear phase.
%
%   Both types of filters are basically lowpass and have a gain of L in
%   the passband.
%
%   Example: Compare two filter with different bandlimitedness factors.
%   b1 = intfilt(5,10,1);  % Signal occupies full Nyquist interval
%   b2 = intfilt(5,10,.8); % Signal occupies 80% of Nyquist interval
%   hfvt = fvtool(b1,1,b2,1); legend(hfvt,'Alpha = 1','Alpha = 0.8');
%
%   See also INTERP, DECIMATE, RESAMPLE.

%   Reference:   Oetken, Parks, and Schuessler, "New Results in the
%      design of digital interpolators," IEEE Trans. Acoust., Speech,
%      Signal Processing, vol. ASSP-23, pp. 301-309, June 1975.

%   Copyright 1992-2013 The MathWorks, Inc.

narginchk(3,3)

% Cast to enforce Precision Rules
R = signal.internal.sigcasttofloat(R,'double','intfilt','L','allownumeric');
L = signal.internal.sigcasttofloat(L,'double','intfilt','P','allownumeric');                                 

if nargin == 3,
    if ischar(freqmult),
        type = freqmult;
        n = L;
    else
        % Cast to enforce Precision Rules
        freqmult = double(freqmult);
        type = 'b';
    end
end 

if strcmp(type(1),'b') || strcmp(type(1),'B'),
        
    n=2*R*L-1;

    if freqmult == 1
        
        M = [R R 0 0];
        F = [0 1/(2*R) 1/(2*R) .5];
    else
        
        M=R*[1 1];
        F=[0 freqmult/2/R];

        for f=(1/R):(1/R):.5,
            F=[F f-(freqmult/2/R) f+(freqmult/2/R)]; %#ok<*AGROW>
            M=[M 0 0];
        end;

        if (F(length(F))>.5),
            F(length(F))=.5;
        end;
    end

    h=firls(n-1,F*2,M);

elseif strcmp(type(1),'l') || strcmp(type(1),'L'),

% Inputs:
%    n   order of polynomial interpolation  (should be ODD!!!!)
%    R   discrete time sampling period (input to filter assumed 0 else)

    if n == 0
        h = ones(1,R);
        return
    end

    if ~(rem(n,2) || rem(R,2))
       warning(message('signal:intfilt:InvalidParamLinPhase', 'R', 'N'));
    end

    t=0:n*R+1;
    l=ones(n+1,length(t));
    for i=1:n+1
        for j=1:n+1
            if (j~=i)
                l(i,:)=l(i,:).*(t/R-j+1)/(i-j);
            end;
        end;
    end;

%    plot(t/R,l')
%    title('Lagrange Polynomials');
%    xlabel('time/R');
%    grid
%    l
%    pause;

    h=zeros(1,(n+1)*R);  

    for i=0:R-1
        for j=0:n
            h(j*R+i+1)=l((n-j)+1,round((n-1)/2*R+i+1));
        end;
    end;

    if h(1) == 0,
        h(1) = [];
    end

else
    error(message('signal:intfilt:InvalidParamFilterType'))
end

 if nargout > 1   
     a = 1;   
 end