Can I "reverse" spectral leakage by solving a linear system?

Question

When I compute a spectrogram of (say) a piece of music, there is a lot of frequency "smearing." Often we can reasonably expect that the "true" generating process is much sparser in frequency (i.e. maybe 10-20 active frequencies instead of hundreds).

We also know that plain old sine waves become smeared in spectrograms.

So, can't we model the corrupting "frequency leakage" at a specific timestep as a linear equation $$ Ax = b $$ Where

• $b$ is the known vector representing the STFFT-computed spectral densities at the various frequencies
• $x$ is the unknown, sparse vector of generating frequencies
• $A$ is the (computable) "smearing matrix" which maps power from sparse generating frequencies to the "smeared" set of frequencies which STTFT produces (We can just compute this by taking the STFFT of various sine wave snippets)

And naturally one could extend this to further encourage sparsity, e.g. $$ \min_{x} \|Ax - b\| + λ \|x\|_1 $$

I have not seen this approach discussed, but I am very new to DSP. Is this a known technique, or is there some reason why it would not work in "sharpening" a spectrogram?

Answer 1

This looks sort of like the standard approach to the maximum likelihood estimator for tones.

Kay (1) starts with the signal $$ x[t] = A \cos(\omega_0 t + \phi) + \epsilon[t] $$ and assumes $\epsilon[t]$ is an independent, identically distributed Gaussian noise source with variance $\sigma_\epsilon^2$.

The likelihood function for this is: \begin{align} \newcommand{\bx}{\mathbf{x}} p(\bx ; \mathbf{\theta}) = \frac{1}{(2\pi \sigma_\epsilon^2)^{T/2} } \exp\Big( -\frac{1}{2\sigma_\epsilon^2} \sum_{t=0}^{T-1} \big(x[t] - \tilde{A} \cos(\tilde{\omega}_0(t-\nu) + \tilde{\phi}) \big)^2 \Big) \end{align} where our parameters of interest are $\mathbf{\theta} = \Big [ \tilde{A}, \tilde{\omega}_0, \tilde{\phi} \Big ]^T$.

Using the log-likelihood, and forming the vector version gives: $$\newcommand{\bc}{\mathbf{c}} \newcommand{\bs}{\mathbf{s}} \newcommand{\bH}{\mathbf{H}} \begin{align} L(\theta') &= \frac{T}{2}\log_e(2\pi \sigma_\epsilon^2) + (\bx -\tilde{\alpha}_c \bc - \tilde{\alpha}_s \bs )^T (\bx -\tilde{\alpha}_c \bc - \tilde{\alpha}_s \bs )\\ &= \frac{T}{2}\log_e(2\pi \sigma_\epsilon^2) + (\bx - \bH\underline{\alpha} )^T (\bx - \bH \underline{\alpha} ) \end{align}$$ where \begin{align} \bx &= \left [ x[0], x[1], \ldots, x[T-1] \right]^T\\ \bc &= \left [ \cos\left(-\tilde{\omega}_0 \frac{T-1}{2}\right), \ldots , \cos(-\frac{\tilde{\omega}_0}{2}), \cos(\frac{\tilde{\omega}_0}{2}), \ldots \cos\left(\tilde{\omega}_0 \frac{T-1}{2}\right) \right]^T\\ \bs &= \left [ \sin\left(-\tilde{\omega}_0 \frac{T-1}{2}\right), \ldots , \sin(-\frac{\tilde{\omega}_0}{2}), \sin(\frac{\tilde{\omega}_0}{2}), \ldots \sin\left(\tilde{\omega}_0 \frac{T-1}{2}\right) \right]^T\\ \bH &= \left [ \bc \ \bs \right ]\\ \underline{\alpha} &= \left[ \tilde{\alpha}_c\ \tilde{\alpha}_s\right ]^T\\ \tilde{\alpha}_c &= \tilde{A} \cos \phi\\ \tilde{\alpha}_s &= -\tilde{A} \sin \phi \end{align} assuming $T$ is even.

Like I said, sort of. The $\underline{\alpha}$ here is analogous to your $x$, but $\underline{\alpha}$ is in no way sparse.

Interesting question. I'll think about it some more.

(1) S. M. Kay, Fundamentals of Statistical Signal Processing: Estimation Theory. Prentice Hall, 1997.