How to obtain Real and Imaginary parts of symbolic polynomial

I have the below code. I cannot figure out how to obtain the real and imaginary parts of the symbolic polynomial (I have already solved the complete problem in standard Matlab, trying to learn the Symbolic Toolbox). My code is:

polynomial = [1 10 45 105 105]; % This polynomial is the Bessel Filter 
% polynomial computed above, and changes with the order of the filter.

p = poly2sym(polynomial, s)

pw = subs(p, 1*j*w)

% attempt from Walter Roberson's comment at: 
% https://www.mathworks.com/matlabcentral/answers/
% 465739-real-and-imaginary-part-of-this-function-please
real_pw = simplify( rewrite(real(pw), 'exp') )
imag_pw = simplify( rewrite(imag(pw), 'exp') )
% NOT WORKING FOR ME

% Assuming I can get the above done, my next steps would be to convert the
% below (from standard Matlab) to symbolic.  Any help with the below would
% also be appreciated.

sq_real = conv(real_pw, real_pw);
sq_imag = conv(imag_pw, imag_pw);
sq_denom = sq_real + sq_imag;

The above code gives me:

real_pw =

10*imag(w^3) - 45*real(w^2) + real(w^4) - 105*imag(w) + 105

imag_pw =

imag(w^4) - 45*imag(w^2) - 10*real(w^3) + 105*real(w)

Instead, I am expecting something like (obtained from my working plain Matlab code):

real_poly = 1 0 -45 0 105

imag_poly = 0 -10 0 105 0

Any help you can provide is appreciated!

Answer 1

[Paul]

From the context of the code, it appear that w is a real number. Assuming it being so will help:

>> syms s
syms w real
polynomial = [1 10 45 105 105]; % This polynomial is the Bessel Filter 
% polynomial computed above, and changes with the order of the filter.
p = poly2sym(polynomial, s);
pw = subs(p, 1j*w);  % note use of 1j, not 1*j
real_pw = simplify( rewrite(real(pw), 'exp') )
imag_pw = simplify( rewrite(imag(pw), 'exp') )

% but this is simpler
real_pw = real(pw)
imag_pw = imag(pw)
 
real_pw =
 
w^4 - 45*w^2 + 105
 
 
imag_pw =
 
- 10*w^3 + 105*w
 
 
real_pw =
 
w^4 - 45*w^2 + 105
 
 
imag_pw =
 
- 10*w^3 + 105*w

I don't believe that conv is supported for symbolic polynomials (version 2019a)

>> conv(real_pw,real_pw)
Error using conv2
Invalid data type. First and second arguments must be numeric or logical.

Error in conv (line 43)
c = conv2(a(:),b(:),shape);
 

But conv is the same as polynomial multiplication, so maybe this is what you want:

>> sq_real = real_pw*real_pw
 
sq_real =
 
(w^4 - 45*w^2 + 105)^2
 
>> sq_imag = imag_pw*imag_pw
 
sq_imag =
 
(- 10*w^3 + 105*w)^2
 
>> sq_denom = sq_real + sq_imag
 
sq_denom =
 
(- 10*w^3 + 105*w)^2 + (w^4 - 45*w^2 + 105)^2
 
>> expand(sq_denom)
 
ans =
 
w^8 + 10*w^6 + 135*w^4 + 1575*w^2 + 11025
 

[Craig]

Very Cool Paul! Thank you very much. And of course, my last step will now be:

w0 = solve(gain^2 / sq_denom - 1/2) % same as solve(gain / sqrt(sq_denom) - sqrt(2)/2)
w0 = w0(w0>0)

I did quite a bit of googling but could not find a way to do the above in only one line ('positive', assume, etc.). If you have a better idea, feel free to comment on it.

[Paul]

sq_denom contains only even powers of w, so you can sub z = w^2 and make it into a fourth order polynomial in z. But that still looks like it will be difficult to solve symbolically. I couldn't get solve to provide much of anything useful, but I'm not that familiar with that function. Maybe someone with more expertise can help. You might need to post as a new question for this part of the problem.

There's a function vpasolve that can find numerical solutions.

Also, sq_denom could be found more directly:

>> simplify(pw*conj(pw))
 
ans =
 
w^8 + 10*w^6 + 135*w^4 + 1575*w^2 + 11025

[Craig]

Thank you very much for the additional comment, clever line of code. I think vpasolve may be the way to go, the equation order goes up with the filter order. This equation is only for a 4 pole filter, 6 and 8 pole filters are not uncommon.

[Walter Roberson]

You do not happen to show the definition for gain so we cannot test this.

Remember that you can add

assume(w >= 0)

and you can also add 'real', true to solve()

[Craig]

Sorry Walter, I should have included:

gain = poly2sym(polynomial(end))

Thank you very much for the advice on 'real', true and on assume(w >= 0). I have incorporated them both into the code and it is working very well now. I appreciate your time and comment!

[Paul]

Walter may be able to work some symbolic magic; he's really good at that. The best I can do is get the numerical answer:

sq_denom(w)=simplify(pw*conj(pw));
gain = polynomial(end);
f(w) = sq_denom(w) - 2*gain^2; % this is the equation to solve for w
fplot(f(w));  % looks like there are two real roots, which makes sense from the form of f(w)
w0 = vpasolve(f,w);
w0(imag(w0)==0)
 
ans =
 
 -2.1139176749042158430196891286339
  2.1139176749042158430196891286339

[Walter Roberson]

solve(f, 'maxdegree', 4) I suspect

If so then the solutions are not likely to clarify anything, as they will be more than 1000 characters each.