Before I state the main problem which I'm interested, I think it's better to provide some background. Suppose that $$w(\xi)=\sum_{n=-N}^Nc_ne^{-in\xi}$$ is a $r\times r$ matrix of trigonometric polynomials, here $c_n$'s are $r\times r$ matrices of real numbers. If $w(\xi)$ is positive semi-definite and $\log(\text{det}(w(\xi))$ is (Lebesgue) integrable on $[0,2\pi]$, then by the operator Fejér-Riesz Theorem, we have the following spectral factorization of $w(\xi)$: $$w(\xi)=\overline{p(\xi)}^Tp(\xi),(*)$$
where $$p(\xi)=\sum_{n=0}^Nd_ne^{-in\xi}$$
is a $r\times r$ matrix of trigonometric polynomials with real coefficients.
Now here comes my question: by using maple, what would be the best way to compute $p(\xi)$ with a given $w(\xi)$?
The most straight forward approach is to parametrize a matrix $p(\xi)$ of trigonometric polynomials, with unknown coefficients $d_n$'s. Then by $(*)$, the coeffifients for each $e^{-in\xi}$ on both sides of the identity have to be the same, and thus we obtain a system of quadratic equations with variables $d_n$'s. Now we can use the "solve" command in maple to directly solve the quadratic system. The major shortcoming of this approach is that, it's highly inefficient, especially when $N$ and $r$ are not that small.
Any suggestion will be greatly appreciated.