0
$\begingroup$

The polynomial here is impulse coefficients of minimum phase FIR filter, or it's impulse response. A code somewhere tries to finding roots of polynomial on the circle. It have roots of polynomial ru and it's derivatives root rdu. The algorithm to finding them uses weird merit to increase unit circle tolerance and that is:

if ~rem(length(ru),2) && length(rdu) == length(ru)/2 && ~isempty(ru)
    No longer increase tolerance and use rdu instead of ru

But why the roots of derivative must be half the original polynomial? Why multiple of two?

Appendix

function [ru,tol,Nsz] = findUnitCircleZeros(b,r,tol,Nsz)
% Find single multiplicity zeros on the unit circle


% Compute the zeros of the derivative
% This will contain the same zeros on the unit circle as b, but not double
rd = roots(polyder(b));

% Initialize flag indicating whether to use the zeros from the derivative
% polynomial or not. Don't use derivative zeros by default.
derivflag = 0;

% Initialize flag
flag = 1;

while tol < 1e-3 && flag
    % Find the zeros on the unit circle
    ru = r(abs(abs(r)-1) < tol);
        
    % Find the zeros on the unit circle corresponding to the derivative
    rdu = rd(abs(abs(rd)-1) < tol);
    
    if nargin > 3
        % Check if number of single zeros found corresponds to number
        % given number of zeros
        if length(rdu) >= Nsz
            rdu = rdu(1:Nsz);
            flag = 0;
            derivflag = 1;
        elseif ~rem(length(ru),2) && length(ru)/2 == Nsz 
            flag = 0;
        else
            % Decrease the tolerance
            tol = 10*tol;
        end
    elseif ~rem(length(ru),2) && length(rdu) == length(ru)/2 && ~isempty(ru)
        % Check that the length of the derivative zeros half of the length of
        % the original zeros on the upper unit circle. If not, decrease
        % tolerance and try again
        flag = 0;
        Nsz = length(rdu);
        if max(abs(polyval(b,ru))) > max(abs(polyval(b,rdu)))
            derivflag = 1;
        end
    else
        % Decrease the tolerance
        tol = 10*tol;
    end
end

% If max tolerance has been reached and number of unit circle zeros has not
% been specified, use the zeros on the unit circle found from the original
% polynomial (not the derivative)
if tol >= 1e-3 && nargin < 4
    Nsz = length(ru)/2;
end  

if ~isempty(ru)
    if derivflag
        % Use unit circle zeros from derivative
        ru = rdu;
    else
        ru = getsinglezeros(ru,Nsz);
    end
    
    if tol < 1e-3
        % Force length of zeros to be exactly one, keep the angle the same
        ru = exp(1i*angle(ru));
    end
else
    Nsz = 0;
end
$\endgroup$
2
  • $\begingroup$ Could you share the complete code? It's a bit hard to say anything useful with such little context. $\endgroup$
    – Matt L.
    Jun 26, 2023 at 13:37
  • $\begingroup$ @MattL. I think it is simple, if you feel any question ask that from me, if I though the deficit is deep then I will add full code. In this way your mind directed exactly toward problem. Though full code added. $\endgroup$ Jun 26, 2023 at 13:39

1 Answer 1

2
$\begingroup$

I'm missing some context but from what I know about such algorithms I can make an educated guess about what's going on. I don't think that the polynomial provided to the function is a minimum-phase polynomial. I believe it is a zero-phase polynomial with a non-negative amplitude response. Since it is non-negative, any zeros on the unit circle must be double (or of a larger even order). Consequently, the number of zeros on the unit circle must be even. Creating a minimum-phase filter from the given response can be done by keeping just one of each double zero on the unit circle. Of course, the zeros that are not on the unit circle must be taken care of somewhere else in the code (keep the ones inside, throw away the ones outside).

The derivative of the polynomial has single roots where the original polynomial has double roots.

For example, if $p(z)$ is

$$p(z) = (z-A)^2(z-B)^2$$

then the derivative is

\begin{align*} \frac{p(z)}{dz} &= 2(z-A)(z-B)^2 + 2(z-A)^2(z-B) \\ &= 2(z-A)(z-B)(z-B + z-A)\\ &= 4(z-A)(z-B)(z - (A+B)/2) \end{align*}

I think that ideally you wouldn't need these roots, but they're used as a second option if there's some (numerical) problem.

$\endgroup$
0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.