Matlab help DSP programs

[BAT_MAN]

Can anybody explain me the logic for the function file although comment are given but i was unable to understand

function Factors = factorize(polyn)
format long; Factors = [];
% Use threshold of 1e-8 instead of 0 to account for
% precision effects
THRESH = 1e-8;
%
proots = roots(polyn); % get the zeroes of the polynomial
len = length(proots); % get the number of zeroes
%
while(len > 1)
if(abs(imag(proots(1))) < THRESH) % if the zero is a real zero
fac = [1 -real(proots(1))];
% construct the factor with proots(1) as zero
Factors = [Factors;[fac 0]];
else % if the zero has imaginary part get all zeroes whose
% imag part is -ve of imaginary part of proots(1)
negimag = imag(proots)+imag(proots(1));
% get all zeroes which have same real part as proot(1)
samereal = real(proots)-real(proots(1));
%find the complex conjugate zero
index = find(abs(negimag) <THRESH & abs(samereal)<THRESH);
if(index) % if the complex conjugate exists
fac = [1 -2*real(proots(1)) abs(proots(1))^2];
%form 2nd order factor
Factors = [Factors;fac];
else % if the complex conjugate does not exist
fac = [1 -proots(1)];
Factors = [Factors;[fac 0]];
end
end
polyn = deconv(polyn,fac);
%deconvolve the 1st/2nd order factor from polyn
proots = roots(polyn); %determine the new zeros
len = length(polyn); %determine the number of zeros
end

 

[just4me]

I think the attachment was missing. Here it is in .rar format

 

[just4me]

hey, It wasn't very difficult to analyze how it works once you have aligned the code properly and then debug it step by step for 1 particular polynomial.

Here's how it works:

1) first we find the roots of the polynomial

2) There is a while loop which runs for each root.

3) If the root is real then obviously the factor is [1 -root] since root of (x-a) is a.

4) If the root is imaginary then the program tries to find another root which is conjugate. If found then it combines these two to form a 2nd degree polynomial.
This is done using the following logic:

(x-(a+ib))(x-(a-ib)) = x^2 - 2ax + (a+ib)(a-ib) = x^2 - 2ax + |a+ib|^2

5) If it does not find the complex conjugate root, then it again declares a factor like [1 (a+ib)] which is the root of (x - (a+ib)).

6) In each iteration, it uses deconv() which finds the remainder of the polynomial when the polynomial is divided by the factor.

Hope this helps in understanding.

I am also attaching the .m file which is aligned properly.