Caucer fraction expansion fraction expansion for filter synthesis

2s^2+2s+1|2s^3 +2s^2 +2s+1 | s

2s^3+2s^2+s

-------------------------------------------------------

s+1|2s^2+2s+1| 2s

So want to know how to pereform Caucer fraction expansion fraction expansion for filter synthesis. For a N/M polynomial. so i get a constant*s mainly not s+1,or s-1 which residue gives sometimes

The example below is the one i am trying to figure out

N = 18889465931478580854784*s^13 + 47283973869184172490752*s^12 + 380144960515421138583552*s^11 + 706526257926373225921184*s^10 + 2291017465069233730748416*s^9 + 2885191461803855148054208*s^8 + 4836216961066520318836736*s^7 + 3887431410705412021850560*s^6 + 3598519003466394276724736*s^5 + 1708028481397562746282496*s^4 + 850144437757157128011776*s^3 + 199942013009514701125133*s^2 + 45117763022779280523264*s + 2734820905916453355520; % Example numerator

M = 47283965291789216579584*s^12 + 117964927714854248120320*s^11 + 706526147746757277255008*s^10 + 1167883982950116513284096*s^9 + 2885191063569152216303936*s^8 + 3033824177148917979283456*s^7 + 3887430981188531521547840*s^6 + 2475385520141079517593600*s^5 + 1708028355427532732161536*s^4 + 587964405289678977105920*s^3 + 199942002958859758602739*s^2 + 26228297618450611699712*s + 2734820905916453355520; % Example denominator want to do it for this

Answer 1

[Star Strider]

The Cauer decomposition is useful. Your post is only the second time I’ve encountered a question about it here. The residue function appears to be the only opttion. I experimented with the Symbolic Math Toolbox, and while is ‘sort of’ works, it only does so with textbook examples.

Try this —

syms s 
H = (2*s^5+40*s^3+128*s)/(s^4+10*s^2+9)
[Hn1,Hd1] = numden(H)
T1 = quorem(Hn1,Hd1)
H2 = Hn1/Hd1 - T1
[Hn2,Hd2] = numden(H2)
T2 =  quorem(Hd2,Hn2)
H3 = Hd2/Hn2 - T2
[Hn3,Hd3] = numden(H3)
T3 = quorem(Hd3,Hn3)
H4 = Hd3/Hn3 - T3
[Hn4,Hd4] = numden(H4)
T4 = quorem(Hd4,Hn4)
H5 = Hd4/Hn4 - T4
[Hn5,Hd5] = numden(H5)
T5 = quorem(Hd5,Hn5)

This gives the same result as in W-K Chen, Passive and Active Filters Wiley 1986 (ISBN 0-471-82352-X), P. 113. (Yes, it’s been a while since I’ve done either Cauer or Foster decompositions.)

I stepped through all this manually, to be sure that it gave the correct result. You can probably automate it with a loop, either using this approach or deconv and residue. The ‘T’ values are the component values.

With a ‘real world’ transfer function, it may not produce such a neat result.

EDIT — There still seems to be something wrong with Answers, in that it does not produce the typeset LaTeX results from the Symbolic Math Toolbox in submitted posts, although they are visible while the code is being run initially. (I have reported this to the powers-that-be in MathWorks.) It doesn’t even produce ASCII results, so you’ll have to run this again here or on your own computer to see the results. This seems to be isolated to Answers, since I get the same result in Firefox and DuckDuckGo, so it seems to be browser-independent.

.

[JLC]

Any sugesstions for an other program to use besides matlab ?

[Star Strider]

No.. I only use MATLAB for these — and every other — problem that MATLAB is appropriate for.

That aside, I suspect all other polynomial continued fraction approaches would give the same result. The Wikipedia artticle on contued fractions (that I would normally expect to proviide guidance if there are any other options available) is curiously silent. It does not specifically discuss polynomial continued fractions.

A specific search for MATLAB applications returned Continued Fractions and Computations however while informative, that deals only with numeric continued fractions, not polynomials. There do not appear to be any File Exchange contributions.

Apparently, what I posted here is about as good as it gets. (I also explored the partfrac function, however while it does a partial fraction expansion useful for inverting Laplace transforms, it is not appropriate for these sorts of problems.)

I cannot figure out from your notation what you want to do, or I would have posted a separate section using your polynoimials, now that I got my original code to work correctly.

EDIT — (4 Oct 2024 at 19:50)

A loop version of this (that can be enclosed a function if desired):

denpoly = [1 1];
syms s 
H = (2*s^5+40*s^3+128*s)/(s^4+10*s^2+9)
k = 0;
H = 1/H;
while numel(denpoly) > 1                            % Loop Version (Can Be A Function)
    k = k + 1;
    % disp("k = "+k)
    [Hn,Hd] = numden(1/H);
    denpoly = sym2poly(Hd);
    Comp(k) = quorem(Hn,Hd);
    H = Hn/Hd - Comp(k);
end

fprintf(['Component Values: \n',repmat('\t\t\t%s\n',1, numel(Comp)), '\n'], Comp)

Component Values: 
			2*s
			s/20
			(40*s)/9
			(9*s)/140
			(70*s)/9

.

[JLC]

N = 18889465931478580854784*s^13 + 47283973869184172490752*s^12 + 380144960515421138583552*s^11 + 706526257926373225921184*s^10 + 2291017465069233730748416*s^9 + 2885191461803855148054208*s^8 + 4836216961066520318836736*s^7 + 3887431410705412021850560*s^6 + 3598519003466394276724736*s^5 + 1708028481397562746282496*s^4 + 850144437757157128011776*s^3 + 199942013009514701125133*s^2 + 45117763022779280523264*s + 2734820905916453355520; % Example numerator

M = 47283965291789216579584*s^12 + 117964927714854248120320*s^11 + 706526147746757277255008*s^10 + 1167883982950116513284096*s^9 + 2885191063569152216303936*s^8 + 3033824177148917979283456*s^7 + 3887430981188531521547840*s^6 + 2475385520141079517593600*s^5 + 1708028355427532732161536*s^4 + 587964405289678977105920*s^3 + 199942002958859758602739*s^2 + 26228297618450611699712*s + 2734820905916453355520; % Example denominator want to do it for this example

[Star Strider]

It is likely not possible to get the correect component values for your transfer fucntion.

This is the best I can do with it —

denpoly = [1 1];
syms s 
% H = (2*s^5+40*s^3+128*s)/(s^4+10*s^2+9)
N =  18889465931478580854784*s^13 + 47283973869184172490752*s^12 + 380144960515421138583552*s^11 + 706526257926373225921184*s^10 + 2291017465069233730748416*s^9 + 2885191461803855148054208*s^8 + 4836216961066520318836736*s^7 + 3887431410705412021850560*s^6 + 3598519003466394276724736*s^5 + 1708028481397562746282496*s^4 + 850144437757157128011776*s^3 + 199942013009514701125133*s^2 + 45117763022779280523264*s + 2734820905916453355520;  % Example numerator
M =  47283965291789216579584*s^12 + 117964927714854248120320*s^11 + 706526147746757277255008*s^10 + 1167883982950116513284096*s^9 + 2885191063569152216303936*s^8 + 3033824177148917979283456*s^7 + 3887430981188531521547840*s^6 + 2475385520141079517593600*s^5 + 1708028355427532732161536*s^4 + 587964405289678977105920*s^3 + 199942002958859758602739*s^2 + 26228297618450611699712*s + 2734820905916453355520;  % Example denominator    want to do it for this example
H = N/M
k = 0;
H = 1/H;
while numel(denpoly) > 1                            % Loop Version (Can Be A Function)
    k = k + 1;
    % disp("k = "+k)
    [Hn,Hd] = numden(1/H);
    denpoly = sym2poly(Hd);
    Comp(k) = quorem(Hn,Hd);
    H = Hn/Hd - Comp(k);
end

fprintf(['Component Values: \n',repmat('\t\t\t%s\n',1, numel(Comp)), '\n'], vpa(Comp,5))

Component Values: 
			0.39948988658017015040968544781208*s + 0.0033452836487053438929706317139789
			0.48496233525588650081772357225418*s + 0.028056395504794551243321620859206
			0.65144734460659492469858378171921*s + 0.15426965412049753467726986855268
			0.55717724103647014999296516180038*s + 0.26010422572608149494044482707977
			0.869940686513359651144128292799*s + 0.47110725868242298020049929618835
			0.74864274422543530818074941635132*s + 0.2593341243966733600245788693428
			1.433139392115435839514248073101*s + 0.46468227734453648736234754323959
			0.94562660091423822450451552867889*s + 0.086912311336121206295501906424761
			2.8071843022230495989788323640823*s + 0.65235246001054747466696426272392
			1.0287572205440937977982684969902*s + 0.02515526541379742297976918052882
			6.3181395311366941314190626144409*s + 0.9203349802396587620023638010025
			1.5125986176108199288137257099152*s + 0.32394380064283723186235874891281
			4.7095560250536436797119677066803*s - 0.88272785194806147046620026230812

The magnitudes of the second part of each vector are not negligible with respect to the magnitudes of the coefficents of s, so it is not possible to simply ignore them. Not all transfer functions can be correctly decomposed into their components.

.

[JLC]

Was thinking about not using qourem and just write it out in code ,like just with multiply/divide and subtract functions.Becasue the functions i am dealing with the numerator is always one order higher .

[Star Strider]

I doubt that you can do this without quorem or its numeeric equivalent (deconv).

Your only other option might be to synthesize and equivalent RL circuit (in thiis instance, having 13 circuit elements) and then seee eif you can manage a least-squarees estimate of its components. (I have neveer tried that.)

syms s 
N =  18889465931478580854784*s^13 + 47283973869184172490752*s^12 + 380144960515421138583552*s^11 + 706526257926373225921184*s^10 + 2291017465069233730748416*s^9 + 2885191461803855148054208*s^8 + 4836216961066520318836736*s^7 + 3887431410705412021850560*s^6 + 3598519003466394276724736*s^5 + 1708028481397562746282496*s^4 + 850144437757157128011776*s^3 + 199942013009514701125133*s^2 + 45117763022779280523264*s + 2734820905916453355520;  % Example numerator
M =  47283965291789216579584*s^12 + 117964927714854248120320*s^11 + 706526147746757277255008*s^10 + 1167883982950116513284096*s^9 + 2885191063569152216303936*s^8 + 3033824177148917979283456*s^7 + 3887430981188531521547840*s^6 + 2475385520141079517593600*s^5 + 1708028355427532732161536*s^4 + 587964405289678977105920*s^3 + 199942002958859758602739*s^2 + 26228297618450611699712*s + 2734820905916453355520;  % Example denominator    want to do it for this example
H = N/M

poles = vpa(solve(M == 0))
zeros = vpa(solve(N == 0))

format long
dpoles = double(poles)

dpoles = 
 -0.079368868777460 + 0.131661886100565i
 -0.079368868777460 - 0.131661886100565i
 -0.138189576136688 + 0.401808290776037i
 -0.138189576136688 - 0.401808290776037i
 -0.170100539607854 - 0.679254719474679i
 -0.170100539607854 + 0.679254719474679i
 -0.215942879707496 + 1.079764865004026i
 -0.215942879707496 - 1.079764865004026i
 -0.310728180925734 + 1.715571246753463i
 -0.310728180925734 - 1.715571246753463i
 -0.333079377752649 - 2.718333920821942i
 -0.333079377752649 + 2.718333920821942i

dzeros = double(zeros)

dzeros = 
 -0.081085573137939 + 0.000000000000000i
 -0.103417324337106 - 0.287851857415961i
 -0.103417324337106 + 0.287851857415961i
 -0.152113684584756 - 0.527943267715307i
 -0.152113684584756 + 0.527943267715307i
 -0.183812958672276 - 0.848877024821345i
 -0.183812958672276 + 0.848877024821345i
 -0.247008399587547 - 1.329699995391562i
 -0.247008399587547 + 1.329699995391562i
 -0.354425990260786 - 2.110936814088112i
 -0.354425990260786 + 2.110936814088112i
 -0.170275222998260 - 3.157834914922368i
 -0.170275222998260 + 3.157834914922368i

figure
scatter(real(dpoles), imag(dpoles), 'x')
hold on
scatter(real(dzeros), imag(dzeros), 'o')
xline(0, '--k')
hold off
grid
axis('padded')
xlabel('Re')
ylabel('Im')
title('Pole-Zero Plot')

There does not appear to be any way to make this work, and produce a simple decomposition. (I did a few experiments to see if that was poissible.)

.

[JLC]

yeah ,thanks thou

[Star Strider]

My pleasure!

I’m not certain what magic to use to perhaps transform or map your transfer function iinto a version that could work. One of the experiments I tried was creating polynomials from the imaginary parts of the poles and zeros, since the Cauer factorisation requires that (no real parts).

This is that experiment, however it takes too long to run here —

syms s 
N =  18889465931478580854784*s^13 + 47283973869184172490752*s^12 + 380144960515421138583552*s^11 + 706526257926373225921184*s^10 + 2291017465069233730748416*s^9 + 2885191461803855148054208*s^8 + 4836216961066520318836736*s^7 + 3887431410705412021850560*s^6 + 3598519003466394276724736*s^5 + 1708028481397562746282496*s^4 + 850144437757157128011776*s^3 + 199942013009514701125133*s^2 + 45117763022779280523264*s + 2734820905916453355520;  % Example numerator
M =  47283965291789216579584*s^12 + 117964927714854248120320*s^11 + 706526147746757277255008*s^10 + 1167883982950116513284096*s^9 + 2885191063569152216303936*s^8 + 3033824177148917979283456*s^7 + 3887430981188531521547840*s^6 + 2475385520141079517593600*s^5 + 1708028355427532732161536*s^4 + 587964405289678977105920*s^3 + 199942002958859758602739*s^2 + 26228297618450611699712*s + 2734820905916453355520;  % Example denominator    want to do it for this example
H = N/M

poles = vpa(solve(M == 0))
zeros = vpa(solve(N == 0))

format long
impoles = poly2sym(double(imag(poles)),s);
imzeros = poly2sym(double(imag(zeros)),s);

H = expand((imzeros * s)) / impoles;
H = vpa(H, 5)

denpoly = [1 1];
% syms s 
% H = (2*s^5+40*s^3+128*s)/(s^4+10*s^2+9)
k = 0;
H = 1/H;
while numel(denpoly) > 1                            % Loop Version (Can Be A Function)
    k = k + 1;
    % disp("k = "+k)
    [Hn,Hd] = numden(1/H);
    denpoly = sym2poly(Hd);
    Comp(k) = quorem(Hn,Hd);
    H = Hn/Hd - Comp(k);
end

fprintf(['Component Values: \n',repmat('\t\t\t%s\n',1, numel(Comp)), '\n'], Comp)

Multiplying the numerator by s is necessary so that the numeerator has a power that is one greater than the denominator polynomial.

.

Answer 2

[Hitesh]

Hi Liam Jacobs,

You can use MATLAB's built-in functions, "quorem" or "residue," to obtain the constant value. Ensure these functions are used within the loop where you continuously update the numerator and denominator of the polynomials. MATLAB's symbolic toolbox can handle high-order polynomials effectively.

Please refer to below code for your reference:

% Define the symbolic variable
syms s;

% Define the numerator and denominator polynomials
% Example: 14th-order numerator with large coefficients
N = 1e5*s^14 + 2e4*s^13 + 3e3*s^12 + 4e2*s^11 + 5e1*s^10 + 6e0*s^9 + ...
    7e5*s^8 + 8e4*s^7 + 9e3*s^6 + 10e2*s^5 + 11e1*s^4 + 12e0*s^3 + ...
    13e5*s^2 + 14e4*s + 15e3;  % Example numerator

M = 2*s^2 + 2*s + 1;  % Example denominator

% Initialize the continued fraction expansion
cauerExpansion = sym([]);

% Perform polynomial division iteratively
while true
    % Perform polynomial division
    [Q, R] = quorem(N, M, s);
    
    % Append the quotient to the continued fraction expansion
    cauerExpansion = [cauerExpansion, Q];
    
    % Check if the remainder is zero
    if R == 0
        break;
    end
    
    % Check if the remainder is a constant
    coeffs = sym2poly(R);
    if length(coeffs) == 1
        fprintf('Remainder is a constant: %s\n', char(R));
        break;
    end
    
    % Update the numerator and denominator for the next iteration
    N = M;
    M = R;
end

% Display the continued fraction expansion
disp('Cauer Continued Fraction Expansion:');
disp(cauerExpansion);

% Format the continued fraction for display
cfString = sprintf('%s', char(cauerExpansion(1)));
for i = 2:length(cauerExpansion)
    cfString = sprintf('%s + 1/(%s)', cfString, char(cauerExpansion(i)));
end
disp('Formatted Continued Fraction:');
disp(cfString);

For more information on "quorem" and "residue" function, kindly refer the following documentation link:

quorem: https://www.mathworks.com/help/symbolic/sym.quorem.html

Residue: https://in.mathworks.com/help/matlab/ref/residue.html

[JLC]

N = 18889465931478580854784*s^13 + 47283973869184172490752*s^12 + 380144960515421138583552*s^11 + 706526257926373225921184*s^10 + 2291017465069233730748416*s^9 + 2885191461803855148054208*s^8 + 4836216961066520318836736*s^7 + 3887431410705412021850560*s^6 + 3598519003466394276724736*s^5 + 1708028481397562746282496*s^4 + 850144437757157128011776*s^3 + 199942013009514701125133*s^2 + 45117763022779280523264*s + 2734820905916453355520; % Example numerator

M = 47283965291789216579584*s^12 + 117964927714854248120320*s^11 + 706526147746757277255008*s^10 + 1167883982950116513284096*s^9 + 2885191063569152216303936*s^8 + 3033824177148917979283456*s^7 + 3887430981188531521547840*s^6 + 2475385520141079517593600*s^5 + 1708028355427532732161536*s^4 + 587964405289678977105920*s^3 + 199942002958859758602739*s^2 + 26228297618450611699712*s + 2734820905916453355520; % Example denominator

want to do it for this ?