TL;DR:
The basis of this response is that matrix valued PSDs imply vector valued quantities that are expressed as Tensors. The Tensor is explored, briefly, to provide a basis to think (and possibly search further, beyond this post, for) practical applications. Practical applications explored here include DWTI-MRI and Quaternion RGBs.
Let $\Phi_1$ and $\Phi_2$ be two matrix-valued power spectral densities.
How does a Power Spectral Density ends up being matrix-valued?
By applying it over a Tensor field or other sequence of Tensors.
What is a Tensor? And why do I need one?
The answer to "What is a Tensor" varies wildly depending on who one dares to ask. That's because different specialities are interested in different aspects and applications of Tensors. From my point of view, I thought I would provide a "necessity" based 2 min definition to cover how do we get to that "matrix valued" part:
- A "scalar field" type of signal, associates a single quantity ($x \in \mathbb{R}^d$ ) with some point in space. That is, independently of the dimensionality ($d$) of the scalar:
- One dimensional signals such as a discrete sampled sound recording ($x[i]$) are represented as a sequence of real ($x \in \mathbb{R}$) samples of sound pressure at distinct points in time. In this case, all we are interested in is how does the Sound Pressure Level vary over time. Notice here that $x$ is assumed to be one continuous recording without "gaps". In this case, successive samples $i, i+1$ are assumed to have been obtained at some $T_s$ time interval called "Sampling Period" (and of course, $i$ is an integer ($i \in \mathbb{Z}$)) . That is, we don't need to convey that $x[i]$ contains successive samples. It is implied.
- Two dimensional signals, such as grayscale images, are still scalars. In visible photography, a sensor (or film surface) records how does light intensity vary within its Field of View. Again, we associate a single measurement ($x[i,j]$) with a point in space ($i,j$).
- Three dimensional signals, such as volumetric data coming out of structural Magnetic Resonance Imaging scanner, are still scalars. In this case, all we are interested in is the concentration of a substance ($x[i,j,k]$) in three dimensional space ($i,j,k$). This concentration implies that $x$ might be more dense at $i_1,j_1,k_1$ but less dense around $i_2,j_2,k_2$ as it happens with fatty tissue around muscles for example.
- Three dimensional signals such as video are still scalars. That is, we associate light intensity $x[i,j,n]$ with a point in 2D space $i,j$ and a point in time $n$.
- Four dimensional signals such as those coming out of functional Magnetic Resonance Imaging are still scalars. In this case, we associate the concentration of a substance $x[i,j,k,n]$ with a point in 3D space $i,j,k$ and a point in time $n$.
But, what do scalar fields look like and why are they not enough?
This is what a three-dimensional scalar field looks like:
And here is a more realistic one. I am showing this one because the effect is more pronounced.
What do we observe? Well, it's not smooth is it? Why is it not smooth? Because it is only composed by little cubes that represent sampled space without any regard for direction. I only know that a particular point in space had some value but I don't know where was this value coming from.
Contrast this previous image with something like this:
This skull looks like it is composed of tiny little surfaces that tell us not only that some point in space $i,j,k$ is occupied by some bone substance but also what the orientation of this substance, in space, is.
To be able to represent this at the point of measurement, we cannot use a scalar field. We need something different. We need something that can encode direction. We need a vector.
Now, someone might say here: Hold on, you can calculate the local slope of a scalar field by obtaining its Gradient. Which, by the way, is a vector field.
And, for some applications, that could be enough but that would also mean that we make a particular assumption about the underlying phenomenon. That is, that "things" flow from high pressures to low pressures. In other words, we would still sample a scalar field of pressure (or density) values and from that infer flow. What if these regions of "high pressure" don't move? (in reality). Or what if only a small amount of their "density" flows?
What if what we need, genuinely, is to be able to not only record the concentration of a substance $x$ at some point in three dimensional space ($i,j,k$) but also its flow towards some direction $\rho, \theta, \alpha$?
Such is the case with Diffusion Tensor Magnetic Resonane Imaging. MRI is a beautiful "sensor" that can return a huge amount of information about a biological sample in its field of view.
It samples the human body directly in the frequency domain and to produce an "image" it transforms the acquired signals, from the frequency domain to the time domain.
Because of the way magnetic fields interact with magnetic objects (even as tiny as atoms), it is possible to selectively "excite" and "listen" for signals at a particular point in three dimensional space $i,j,k$.
This is how structural and functional MRI works which return scalar fields. You "excite" a particular region and then listen to the signal of atoms (from that particular region) "reacting" to the excitation (typically in the region of MHz, e.g. for H the frequency is ~43 MHz per Tesla). The more atoms you "listen" the higher the $x[i,j,k]$ at that point. Or more generally, $x$ is proportional (even if negatively) to the concentration of specific atoms at that $i,j,k$
But, by modifying the way atoms are excited you can also infer if and how did they move. How do you do that? You "excite" a particular region but "listen" from another region in quick succession. In other words, you "tag" some atoms, wait and then look where did the tag go (but very quickly because the "tag" is "wearing out" with time).
The geometrical orientation of these "excite" / "listen" regions now forms a vector. If you found "a lot of tagged atoms" towards a particular direction then this direction has higher magnitude of "flow", otherwise it does not.
So, in this case, we need to represent both the magnitude of the quantity ("How much substance moved?") and its direction ("In what direction did it move?").
We do this using a Tensor Field. It looks like this:
$$A[i,j,k] = \begin{bmatrix}
A_{ii} & A_{ij} & A_{ik} \\
A_{ji} & A_{jj} & A_{jk} \\
A_{ki} & A_{kj} & A_{kk} \\
\end{bmatrix}
$$
So, imagine a three dimensional grid such as this one:
This grid has a "node" (a point where three axes meet) at distinct locations $i,j,k$. But this time around, instead of having a single number associated with that $i,j,k$ we have a Tensor ($A$).
The MRI populates this matrix by obtaining six sets of measurements corresponding to "geometrical displacements" (i.e. directions) $A_{ij}, A_{ik}, A_{ji}, A_{jk}, A_{ki}, A_{kj}$. Pairs of $indexVar_1 indexVar_2$ denote directions. So $ij$ can be left to right, $ji$ right to left, $ik$ top to bottom, $ki$ bottom to top, $jk$ front to back $kj$ back to front WITH RESPECT TO $[i,j,k]$.
It is important to understand that $A$ is a matrix and it is one thing. Just as some vector is a "packet" (tuple) of two numbers ($\vec{a} = (b,c))$ but we treat it effectively as one thing, in exactly the same way, $A$ is one thing. One, indivisible, inseparable object. ONE thing.
This is important for what is about to follow.
Now, if you check back our $A$, you will notice that it contains a huge amount of information. So much in fact, that we can infer flow towards any direction.ANY direction. No, seriously, look back at $A$ again. $A$ (and of course, not $A$ but the MRI itself) captures "symmetric" flows. It captures flow from left to right and at the same time it captures a flow from right to left. Therefore, we can "see" exchange phenomena. Something travelling from left to right (perhaps mostly) but at the same time another flow (perhaps a low amount) travelling backwards in the same direction at the same site.
Most of the times however what we are interested in is just one vector out of all the possible directions that $A$ can encode (The vector of "most flow"). We extract this via the eigenvectors and eigenvalues of $A$.
The last step in this process is to visualise the flow. In this case, we simply "follow the arrows" to reconstruct an orbit. This orbit coincides with neuronal axons (in the brain) and this is what the DTI images depict. About this image, it is important to understand that this is a static image. This is a snapshot of the tracts, although it has been created by an extremely large amount of measurements that may have taken some time. If we did DTI over time, then we would have one $A$ associated with each $i,j,k,n$.
OK, that's fine, but how do you end up with PSDs that look like matrices?
Remember that $A$ that is one thing? Just like some $x[i]$ (e.g. some integer number) is one inseparable thing?
Find the periodogram of $A[i,j,k]$.
This is where it pays off to think about $A$ as one thing. Because you don't treat it differently than any other vector, although, yes, it certainly is not "just a vector".
Now, I will use Welch's method of deriving the periodogram here which specifies splitting the "signal" into overlapping segments, taking the DFT of those segments, deriving the spectrum and then averaging those segments to derive the final power spectrum.
BUT!!!
The most challenging thing here is NOT the computational part. The most challenging thing here is the connection with reality. What does it mean? What does it look like? In fact, the most challenging thing in the aforementioned process of deriving the periodogram is not the "indexing" (chopping the Tensor field down to overlapping segments) because that's just like slicing a scalar array. Only that now, each one of the entries of this Tensor field is a Tensor. So, this would look like an Array of Tensors, which we cut down into overlapping segments. It is not the averaging, because, averaging is averaging. Tensor addition works as you would expect, you just add the individual components.
The most challenging part, in terms of comprehension is The Discrete Fourier Transform. Again, it's not about calculating it, because, the Discrete Fourier Transform of a Tensor field is the Discrete Fourier Transform of the individual components of the Tensors that compose it. In terms of notation, you still have sinusoids but they are packed in a Tensor "package" this time.
The challenging bit here is the meaning of it and what does it look like given that it describes so much information.
One (useful) view would be to treat it like the DFT of complex data. Because in that case, it's not $x \in \mathbb{R}$ anymore but rather $x \in \mathbb{C}$ and therefore $x$ is now actually a vector (i.e., it has magnitude and direction).
When $x \in \mathbb{C}$ then what the Fourier transform actually does is decompose $x$ to the sum of a set of "squiggly" Lissajous symbols. That is, the complete repertoire of "doodles" that two interacting (orthogonal) sinusoids can create, including a perfect circle (at DC) but more importantly ellipsoids at various orientations.
This concept can be applied as is to the Fourier Transform of the Eigenvector of $A$ described above. In that case what we would be getting after $\mathcal{F}(\lambda \cdot v), A = \lambda \cdot v$ would be a decomposition of $v$ down to a set of "squiggly" three dimensional surfaces including a perfect sphere and also ellipsoids with various orientations (and everything in between).
I suppose that something similar can be obtained to interpret $\mathcal{F}(A)$ although it probably does not look exactly the same because each frequency of the PSD is now a complex matrix and returns a huge amount of directional information about that particular point (which we kind of "discard" when we focus only on the eigenvector).
...a practical application in which the considered spectra are matrix-valued functions
This has been partially covered by the DTI and discussion in the previous section. The only thing that I would like to add is that in terms of interpolating, you can interpolate incredibly complex curves in surprising ways with results similar to this one.
And also that in the signal / image processing space, there is a "Tensor" expression of the DFT that would return spectra you can work with but I am not really sure about interpolation.
Perhaps you can define some sort of "smoothing" filter. That is, given a Tensor field with missing spectral values, fill them in by interpolating the known values around it.
Other than this, I suppose that you can define a metric of "distance" between two spectra (?) to "match" tensor fields.
Beyond these applications in DSP, the next stop would be the Physics SE. They use Tensors a lot in cosmology. I don't quite understand this stuff,...I mean, what does it look like...and sometimes the "why". I guess they like tensors because they express that oneness required to describe space and time...but how do you go from that to, "maybe "time" does not exist", is impenetrable for me.
Anyway, this is why I asked "what is it that you do"? What is your "necessity" for ending up with Matrix valued PSDs? What is the "problem"?
Hope this helps.
