I know we can do (Discrete) Fourier Transform of vectors of polynomial coefficients. This is useful for example when multiplying polynomials, since convolution turns into multiplication in Fourier domain which reduces computational complexity.
But, what if we have a matrix where each element is a polynomial?
Is there some useful sense we can take the Fourier transform of this matrix?
Your question is unclear. With $F$ the DFT and $g\in \Bbb{C}[t]_{n-1}$ (the polynomials of degree $\le n-1$) then $F[g]_m = g(e^{-2i \pi m/n})$.
If $g,h,gh \in \Bbb{C}[t]_{n-1}$ then $F[gh]_m= F[g]_mF[h]_m$.
You are asking if something similar exists for $G \in M_k(\Bbb{C}[t]_{n-1})$, then yes it is $$F[G]_m = G(e^{-2i\pi m/n}) \in M_k(\Bbb{C})$$
With $\times$ the multiplication of matrix, if $G,H,G\times H \in M_k(\Bbb{C}[t]_{n-1})$ then $$F[G \times H]_m = F[G]_m \times F[H]_m$$
We still have the inverse DFT $$G(t) = \sum_{l=0}^{n-1} t^l\frac1n \sum_{m=0}^{n-1} G(e^{-2i\pi m/n})e^{2i \pi lm/n} $$
If you want to unify all the different degrees you'll need to look at $F[G](x) = G(e^{-2i \pi x}), G(t) = \sum_{l=0}^\infty t^l \int_0^1 G(e^{-2i\pi x}) e^{2i \pi l x}dx$