Continue to Site

Welcome to EDAboard.com

Welcome to our site! EDAboard.com is an international Electronics Discussion Forum focused on EDA software, circuits, schematics, books, theory, papers, asic, pld, 8051, DSP, Network, RF, Analog Design, PCB, Service Manuals... and a whole lot more! To participate you need to register. Registration is free. Click here to register now.

Register Log in

factorizing polynomials into as^2 + bs + c for chebyshev filters Matlab

Status
Not open for further replies.
O

oshaye3

Member level 3
Joined
Aug 4, 2010
Messages
62
Helped
0
Reputation
0
Reaction score
0
Trophy points
1,286
Location
london
Activity points
1,822
Hi all,

Thanks to you guys in this forum. You are doing excellent things! I wanted to break my high order polynomials into second orders using Matlab in the form of as^2 + bs + c.

for example : I have a transfer function as
1
----------------------------------------------
s^4 + 1.804 s^3 + 2.627 s^2 + 2.026 s + 0.8286

I need to factorize it into second order form as
1
-----------------------------------
(as^2 + bs + c)(as^2 + bs + c)

Any response is highly appreciated!
 

If you develop (as^2 + bs + c)(as^2 + bs + c) you will obtain:

a^2s^4 +2abs^3 +(2ac + b^2)s^2 + 2bcs + c^2

Then equating to your coefficients

a^2 = 1
2ab = 1.804
2ac+b^2 = 2.627
2bc = 2.026
c^2 = 0.8286

If all the coefficients are positive, then

a=1
b=1.804/2 = 0.902
c=sqrt(0.8286) = 0.91

using these values 2ac+b^2 = 2.63
but 2bc = 1.64 instead of 2.026

so you can't exactly factorize the original 4th order polynomial in a square of a 2nd order polynomial. Do you want, instead, to calculate a more or less rough approximation of the original polynomial ?
 

No, I want to factorise it into 2nd order. There should be a command for that. I am still searching through.
 

With MatLab:
r= roots([1,1.804,2.627,2.026,0.8286 ]);
p1 = poly([r(1),r(2)])
p2 = poly([r(3),r(4)])
​

Without MatLab:
Write the poly s^4+1.804 s^3+2.627 s^2+2.026 s+0.8286 in WolframAlpha and read results.

(s^2+0.528251 s+1.33015) (s^2+1.27575 s+0.62294)​
 

Eduardo, it is very good and helpful alot. What about if someone have unspecified number of order of polynomials. Is there anyway to manipulate the code to factorize the polynomials in Matlab? Thanks
 

An example of m-file could be:

Code: 
function q = polyfact(A)
r = roots(A);
l=1;k=1;
while k <= length(r)
    if isreal(r(k))
        q(l,:) = [0,1,-r(k)]; 
    else 
        q(l,:) = [1,-2*real(r(k)),r(k)*r(k+1)];   
        k=k+1;
    end 
    l=l+1;k=k+1;
end
return

Then, if the poly has some real roots, the output will be:

>> A= [1,1,1,1,1,1,1,-1.804,2.627,2.026,-0.8286 ];
>> polyfact(A)

ans =

1.0000 2.3259 1.6596
1.0000 0.8518 1.7000
0 1.0000 0.7341
1.0000 -0.9692 1.4005
1.0000 -1.6360 0.9320
0 1.0000 -0.3065


that mean
s^2 + 2.3259 s + 1.6596
s^2 + 0.8518 s + 1.7000
s + 0.7341
s^2 - 0.9692 s + 1.4005
s^2 - 1.6360 s + 0.9320
s - 0.3065
 
Thanks Eduardo, It works as a function. Another question is that if I want to factorize them as transfer function such as the code here, it is not given me the results. Look at the code below and the results
code:

H1 = tf(1,den2) % Transfer function
A=den2
den3=polyfact(A)
H2=tf(1,1:den3)

Results:

Transfer function:
1
----------------------------------------------
s^4 + 1.804 s^3 + 2.627 s^2 + 2.026 s + 0.8286


A =

1.0000 1.8039 2.6270 2.0257 0.8286


den3 =

1.0000 0.5283 1.3301
1.0000 1.2755 0.6230

Transfer function:
1
 

oshaye3 said:
...Another question is that if I want to factorize them as transfer function such as the code here, it is not given me the results.
Click to expand...

To do that, polyfact must return a cell array.

Code: 
function q = polyfact(A)
r = roots(A);
l=1;k=1;
while k <= length(r)
    if isreal(r(k))
        q(l) = {[1,-r(k)]};
    else
        q(l) = {[1,-2*real(r(k)),r(k)*r(k+1)]};
        k=k+1;
    end
    l=l+1;k=k+1;
end
return

And the numerator must be another cell array of 1's.

Code: 
den2= [1,1,1,1,1,1,1,-1.804,2.627,2.026,-0.8286 ];
H1 = tf(1,den2)  

p= polyfact(den2);
num = num2cell(ones(size(p)));
H2= tf(num,p)
 
Thanks, Eduardo, I tried it but with the error! check the code and results

code:

H1 = tf(1,den2) % Transfer function
A=den2
p=polyfact(A)
num = num2cell(ones(size(p)));

Results:

A =

1.0000 1.7441 2.7709 2.3972 1.4357 0.4096


p =

1.0000 0.3331 1.1950
1.0000 0.8720 0.6360
0 1.0000 0.5389

??? Error using ==> checkNumDenData at 19
Numerators and denominators must be specified as non-empty row vectors.

Error in ==> testing at 69
H2= tf(num,p)
 

The polyfact function has changed --> Update it (msg #8 )
 
Thanks very much, I have changed the function. it is working,this is the result!
code:

H1 = tf(1,den2) % Transfer function
A=den2
p=polyfact(A)
p_size=ones(size(p))
num = num2cell(p_size);
H2= tf(num,p)

result:

Transfer function:
1
----------------------------------------------------------
s^5 + 1.744 s^4 + 2.771 s^3 + 2.397 s^2 + 1.436 s + 0.4096


A =

1.0000 1.7441 2.7709 2.3972 1.4357 0.4096


p =

[1x3 double] [1x3 double] [1x2 double]


p_size =

1 1 1


Transfer function from input 1 to output:
1
----------------------
s^2 + 0.3331 s + 1.195

Transfer function from input 2 to output:
1
---------------------
s^2 + 0.872 s + 0.636

Transfer function from input 3 to output:
1
----------
s + 0.5389
 

I have ploted the real and Imaginary Poles, the Impulse and step response. I can not plot the frequency response in s domain using freqs



ww = logspace(-1,1);
freqs(H1,ww)
It gives error : Error in ==> testing at 95
freqs(H1,ww)
 

oshaye3 said:
I have ploted the real and Imaginary Poles, the Impulse and step response. I can not plot the frequency response in s domain using freqs
Click to expand...

freqs does not accept cell arrays, only row vectors.

ww = logspace(-1,1);
freqs(num{1},H1{1},ww) % index 1,2,3...
 
It is still giving error! Any suggestion. The error below:

??? Error using ==> lti.subsref at 59
Cell contents reference from a non-cell array object.

Error in ==> testing at 95
freqs(num{1},H1{1},ww) % index 1,2,3...
 

Sorry, I assumed what that H1 was the value returned by polyfact (a cell array)

But if your code was:

Code: 
den2= [1 1.744 2.771 2.397  1.436  4096];
H1 = tf(1,den2) % Transfer function
A=den2
p=polyfact(A)
p_size=ones(size(p))
num = num2cell(p_size);
H2= tf(num,p)
Then H1 and H2 are TF objects! You cannot transfer arguments in any format, freqs accept row vectors only.

To use freqs, you have not need create a TF object (H1,H2)

Code: 
ww = logspace(-1,1);
freqs(1,den2,ww)      % for H1
freqs(num{1},p{1},ww) % for H2  , index 1,2,3...
 
Thanks Eduardo, I have been able to plot the Magnitude and Phase. I have converted the rad/s to Hz by ww/2* pi. I am not sure if the plot is right. I attached the plot with this message. This is the code for the plot

ww = logspace(-1,1);
f_ww= ww/(2*pi); % Convert to Hertz
h_norm = freqs(1,den2,f_ww);
mag = abs(h_norm);
mag = 20*log10(mag);
phase = angle(h_norm);


subplot(2,1,1), loglog(ww,mag) %figure % for H1

This is for 6th order attached !
Another question is that if frequecy is giving, how can we incorporate it to scale the graph eg w=2*pi * f where f is the cutt off frequency. Can we scale it with cutt off freq (f) as well. any tip. Thanks!
 

oshaye3 said:
...This is the code for the plot

ww = logspace(-1,1);
f_ww= ww/(2*pi); % Convert to Hertz
h_norm = freqs(1,den2,f_ww);
mag = abs(h_norm);
mag = 20*log10(mag);
phase = angle(h_norm);
subplot(2,1,1), loglog(ww,mag) %figure % for H1

This is for 6th order attached !
Click to expand...
:? This code not generate that plot.

ww = logspace(-1,1);
f_ww= ww/(2*pi); % Convert to Hertz
h_norm = freqs(1,den2,f_ww);
mag = abs(h_norm);
% mag = 20*log10(mag); <<-- Why? If you plot whith loglog
phase = unwrap(angle(h_norm)); % <<-- type 'help unwrap' in MatLab
subplot(2,1,1), loglog(ww,mag) %figure % for H1

Another question is that if frequecy is giving, how can we incorporate it to scale the graph eg w=2*pi * f where f is the cutt off frequency. Can we scale it with cutt off freq (f) as well. any tip.
Click to expand...


[num,den]= cheby1(6,.5,1,'s'); % Wc = 1 [rad/s]
ww = logspace(-1,1);
h_norm = freqs(num,den,ww);
mag = abs(h_norm); mag = 20*log10(mag);
phase = unwrap(angle(h_norm))*180/pi;

% to display axis with fc [Hz]
fc = 1500 ;
f_ww= ww*fc ;
subplot(2,1,1), semilogx(f_ww,mag)
subplot(2,1,2), semilogx(f_ww,phase)
 
Thanks for your time to guide me and writing this program.

I am trying to figure out the calculation of unknown components in the Multiple feedback. if the transfer function is H(s)=
1
---------------------------------------------------- Let us use only second order sections of transfer function as
R1R3C2C5s^2 + s(R1C5 + R3C5 + R1R3C5/R4) +R1/R4

1
---------------------
s^2 + 0.872 s + 0.636

C_A=1 : C_A=R1R3C2C5
C_B= 0.872 : C_B=(R1C5 + R3C5 + R1R3C5/R4)
C_C=0.636 : C_C=R1/R4

I am using an active topology the link is Electronic filter topology - Wikipedia, the free encyclopedia

How can I write the program to calculate the value of the 3 resistors and two capacitors. Is it going to be writing functions to get the components values?

Let say we have a 3 functions now. Please any tip or is not possible.

C_A= f1( R1 R3 C2 C5)
C_B= f2( R1 R3 R4 C5)
C_C= f3( R1 R4)

Thanks,
 

Status
Not open for further replies.

Similar threads

Part and Inventory Search

Welcome to EDABoard.com

Sponsor

Back
Top