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

    function [nfilt,ff,grid,des,wt,ftype,sign_val,hilbert,neg] = firpminit(order, ff, aa, varargin)
%FIRPMINIT   

%   Copyright 2006-2013 The MathWorks, Inc.

narginchk(3,6);

if all( ff(2:2:end)-ff(1:2:end) == 0),
	error(message('signal:firpminit:InvalidFreqVecZeroWidth'))
end

if order < 3
    error(message('signal:firpminit:InvalidRangeOrder'));
end
%
% Define default values for input arguments:
%
ftype = 'f';
wtx = ones(fix((1+length(ff))/2),1);
lgrid = 16;   % Grid density (should be at least 16)
%
% parse inputs and alter defaults
%
%  First find cell array and remove it if present
for i=1:length(varargin)
    if iscell(varargin{i})
        lgrid = varargin{i}{:};
        if lgrid < 16,
            warning(message('signal:firpminit:InvalidParamGridDensity'));
        end
        if lgrid < 1,
            error(message('signal:firpminit:MustBePositive'));
        end
        varargin(i) = [];
        break
    end
end
if length(varargin) == 1
    if ischar(varargin{1})
        ftype = varargin{1};
    else
        wtx = varargin{1};
    end
elseif length(varargin)==2
    wtx = varargin{1};
    ftype = varargin{2};
end

if isempty(ftype), ftype = 'f'; end

%
% Error checking
%
if rem(length(ff),2)
    error(message('signal:firpminit:MustBeEven'));
end
if any((ff < 0) | (ff > 1))
    error(message('signal:firpminit:InvalidRangeFreqs'));
end
df = diff(ff);
if (any(df < 0))
    error(message('signal:firpminit:InvalidFreqVecNonDecreasingFreqs'));
end
if length(wtx) ~= fix((1+length(ff))/2)
    error(message('signal:firpminit:InvalidDimensions'));
end

if (any(sign(wtx) == 1) && any(sign(wtx) == -1)) || any(sign(wtx) == 0),
    error(message('signal:firpminit:InvalidRangeWeights'));
end

% Determine "Frequency Response Function" (frf)

% The following comment, MATLAB compiler pragma, is necessary to allow the
% Compiler to find the FIRPMFRF private function. Don't remove.
%#function firpmfrf
%#function multiband
%#function lowpass
%#function highpass
%#function bandpass
%#function bandstop
%#function invsinc
%#function hilbfilt
%#function differentiator

if ischar(aa)
    frf = str2func(aa);
    frf_params = {};
elseif isa(aa, 'function_handle');
    frf = aa;
    frf_params = {};
elseif iscell(aa)
    frf = aa{1};
    if ischar(frf)
        frf = str2func(frf);
    end

    frf_params = aa(2:end);
else
    % Check for valid filter length. Ideally we would check in all cases,
    % for now we only do it when aa is a vector
    
    % Cast to enforce precision rules
    aa = double(aa);
    
    exception = 0;
    if any(strncmpi(ftype, {'differentiator','hilbert'}, length(ftype))),
        exception = 1;
    end
    [order,msg1,msg2,msgobj] = firchk(order,ff(end),aa,exception);
    if ~isempty(msg1)
      error(msgobj); 
    end
    
    if ~isempty(msg2),
        warning(message('signal:firpminit:OrderIncreased', ...
          msg2,'h','firpm(N,F,A,W,''h'')'));
    end
    % Grid and weights will be cast to double in firpm to enforce precision
    % rules.
    frf = @firpmfrf;
    frf_params = { aa, strcmpi(ftype(1),'d') };
end

% We need the following code to generate fcn handles to private functions.
% Otherwise the syntax h=firpm(n,f,{@firpmfrf,m},w); will not work
fcns = functions(frf);
if isempty(fcns.file)
    frf = str2func(func2str(frf));
end

%
% Determine symmetry of filter
%
sign_val = 1.0;
nfilt = order + 1;        % filter length
nodd = rem(nfilt,2);      % nodd == 1 ==> filter length is odd
% nodd == 0 ==> filter length is even

hilbert = 0;
if ftype(1) == 'h' || ftype(1) == 'H'
    ftype = 3;  % Hilbert transformer
    hilbert = 1;
    if ~nodd
        ftype = 4;
    end
elseif ftype(1) == 'd' || ftype(1) == 'D'
    ftype = 4;  % Differentiator
    sign_val = -1;
    if nodd
        ftype = 3;
    end
else
    % If symmetry was not specified, call the fresp function
    % with 'defaults' string and a cell-array of the actual
    % function call arguments to query the default value.
    try
        h_sym = frf('defaults', {order, ff, [], wtx, frf_params{:}} );
    catch
        h_sym = 'even';
    end

    if ~any(strcmp(h_sym,{'even' 'odd'})),
        error(message('signal:firpminit:InvalidParamSym', h_sym, frf, '''even''', '''odd'''));
    end

    switch h_sym
        case 'even'
            ftype = 1;  % Regular filter
            if ~nodd
                ftype = 2;
            end
        case 'odd'
            ftype = 3;  % Odd (antisymmetric) filter
            if ~nodd
                ftype = 4;
            end
    end
end

if (ftype == 3 || ftype == 4)
    neg = 1;  % neg == 1 ==> antisymmetric imp resp,
else
    neg = 0;  % neg == 0 ==> symmetric imp resp
end


%
% Create grid of frequencies on which to perform firpm exchange iteration
%
grid = firpmgrid(nfilt,lgrid,ff,neg,nodd);
while length(grid) <= nfilt,
    lgrid = lgrid*4;  % need more grid points
    grid = firpmgrid(nfilt,lgrid,ff,neg,nodd);
end
%
% Get desired frequency characteristics at the frequency points
% in the specified frequency band intervals.
%
% NOTE! The frf needs to see normalized frequencies in the range
% [0,1].
[des,wt] = frf(order, ff, grid, wtx, frf_params{:}); 

%--------------------------------------------------------------------------
function grid = firpmgrid(nfilt,lgrid,ff,neg,nodd)
% firpmgrid
%    Generate frequency grid
nfcns = fix(nfilt/2);
if nodd == 1 && neg == 0
    nfcns = nfcns + 1;
end
grid(1) = ff(1);
delf = 1/(lgrid*nfcns);
% If value at frequency 0 is constrained, make sure first grid point
% is not too small:
if neg ~= 0 && grid(1) < delf
    % Handle narrow bands
    if ff(1) > sqrt(eps)
       grid(1) = ff(1);
    elseif delf < ff(2)
        grid(1) = delf;
    else
        grid(1) = 0.5*(ff(2) - ff(1));
    end
end
j = 1;
l = 1;
while (l+1) <= length(ff)
    fup = ff(l+1);
    newgrid = (grid(j)+delf):delf:(fup+delf);
    if length(newgrid) < 11
        delf1 = ((fup+delf) - (grid(j)+delf)) / 10;
        newgrid = (grid(j)+delf1):delf1:(fup+delf1);
    end
    grid = [grid newgrid];
    jend = length(grid);
    if jend > 1,
        grid(jend-1) = fup;
        j = jend;
    else
        j = jend + 1;
    end

    l = l + 2;
    if (l+1) <= length(ff)
        grid(j) = ff(l);
    end
end
ngrid = j - 1;
% If value at frequency 1 is constrained, remove that grid point:
if neg == nodd && (grid(ngrid) > 1-delf)
    if ff(end-1) < 1-delf
        ngrid = ngrid - 1;
    else
        grid(ngrid) = ff(end-1);
    end
end
grid = grid(1:ngrid);



% [EOF]