What you're asking for is the application to quantum theory of time-frequency analysis (along with time-frequency representations), which are common staples in Digital Signal Processing that haven't quite yet found their way into theoretical physics, with a few notable exceptions, such as: Kaiser and G.H. Battle. Usually, time-frequency analysis is applied to time series - hence the name - but it could also be applied to signal domains of two or more dimensions, like the two-dimensional (or even three-dimensional) signals in graphics and image processing. It can also be applied to developing state-space representations in quantum theory. I talked with Kaiser briefly a while back and he says he's working on a sequel to the work in his research, including the 2011 book that I linked to above, that's more geared toward quantum applications.
The Wigner Function is an example of a time-frequency distribution, but that's really only the tip of the iceberg. There are a lot more time-frequency distributions and transforms than that, and Quantum Physicists are significantly behind the curve on this matter and, as a rule, the theoretical Physics community has Dunning-Kruger on matters related to Fourier Analysis and its extension to Time-Frequency and Time-Scale Analysis. (To put it more directly: when they speak with authority on complementarity, they're outside their field and need to get back to their lane.)
More appropriate in place of time-frequency analysis, actually, is what would be considered "Time Scale" analysis - where in place of frequency on a linear scale, is the version of it on a logarithmic scale, which linked to what's called "scale". This is where octaves are equally-spaced. It is also more closely aligned to how rescaling is approached in renormalization theory. The basis functions in time-scale representations are called "wavelets", though the term is usually reserved for the special case of time-series transforms where the inverse transform has the same form as the forward transform. These are called the Wavelet Transforms - which can be done discretely or continuously. It has heavy application in multi-dimensional signal analysis, specifically image-processing and motion-processing. Battle and Kaiser work primarily with time-scale analysis, rather than time-frequency analysis, in the application of their works to quantum theory.
The Wigner function is the time-frequency transform that results if you apply the operator correspondence to the delta function, itself:
$$δ\left(q - q_0, p - p_0\right) → δ\left(\hat{q} - q_0, \hat{p} - p_0\right),$$
the Wigner transform, itself, being the result of applying this to the wave function representation of a state,
$$W_ψ\left(q_0, p_0\right) = ❬ψ|δ\left(\hat{q} - q_0, \hat{p} - p_0\right)|ψ❭.$$
The main advantage of the Wigner function is that it gives you the obvious results for reasons that are made clear by this way of defining it, e.g.
$$\begin{align}
\int W_ψ\left(q_0, p_0\right) F\left(q_0\right) dq_0 dp_0
&= \int ❬ψ|δ\left(\hat{q} - q_0, \hat{p} - p_0\right)|ψ❭ F\left(q_0\right) dq_0 dp_0\\
&= ❬ψ|\left(\int δ\left(\hat{q} - q_0, \hat{p} - p_0\right) F\left(q_0\right) dq_0 dp_0\right)|ψ❭\\
&= ❬ψ|F(\hat{q})|ψ❭.
\end{align}.$$
Similarly,
$$\int W_ψ\left(q_0, p_0\right) G\left(p_0\right) dq_0 dp_0 = ❬ψ|G(\hat{p})|ψ❭.$$
A "quantized delta function" is ambiguous in that different operator ordering conventions may lead to different implementations of it. In fact, the quantized delta function could, itself, be considered as the very kernel of the operator-orderer that defines the operator ordering convention. The one that goes with the Wigner function is the "Weyl ordering", which symmetrizes the $p$'s and $q$'s, e.g. $\widehat{(qp)} = ½\left(\hat{q}\hat{p} + \hat{p}\hat{q}\right)$. Other variants are those corresponding to other conventions, such as putting all the $p$'s on one side, all the $q$'s on the other, which are instances of the time-frequency distributions known as the Rihaczek distributions (an example of an application).
They each have the disadvantage of being bi-linear, not linear. The more accurate way to write Wigner, for instance, would be as
$$W_{ψψ'}\left(q_0, p_0\right) = ❬ψ|δ\left(\hat{q} - q_0, \hat{p} - p_0\right)|ψ'❭.$$
It's quadratic in $ψ$, with $W_ψ$ actually being $W_{ψψ}$. So, under addition, you get cross-terms:
$$W_{ψ+ψ'} = W_ψ + W_ψ' + W_{ψψ'} + W_{ψ'ψ}.$$
Therefore, Wigner distributions are laden with interference.
Another disadvantage is that the distribution can be negative. In effect, it is an "over-sharpened" picture: it's what you get if you try to de-blur an already sharp image. In fact, if you smear a Wigner distribution with a Gaussian that has a one-SD ellipse sufficiently large, the result will be a strictly non-negative distribution. In turn, the result will be equivalent to a representation by "coherent states".
With either time-frequency or time-scale distributions it is possible to get most of the advantages of Wigner, without the interference or negativity issue; with the transforms being linear, not bi-linear. Wavelets already get you part of the way there, but you can actually go further by just simply dropping the requirement that the inverse transform have the same form as the forward transform.
Example
Here's an example of a Wigner-like time-scale transform constructed in my secret subterranean dungeon Scalographic demo. The mathematics for it are in the description, but I'll replicate the analysis here, to serve as a point of reference.
In the description, the functions $f(t)$ play the role of your wave function $ψ(t)$, while $ψ$ and its Fourier transform $Ψ$ are used to denote the wavelet functions. Sorry for any confusion of notation.
If $f(t)$ is the function for the waveform,
$$f(t,p) = \int f\left(t + \frac{λ}{p}\right) ψ(λ)^* dλ$$
describes its scalogram, with $t$ as the horizontal coordinate, $p$ as the vertical coordinate (displayed on a logarithmic scale), and where $(⋯)^*$ denotes complex conjugate. The original waveform is recovered as
$$f(t) = \int f(t,p) d \log p,$$
which requires that $\int Ψ(γ) d \log γ = 1$ (or more generally, that it be finite and non-zero), where $Ψ(γ) = ∫ ψ(λ) \exp(-2πiγλ) dλ$ is the Fourier transform of $ψ(λ)$.
The scalogram can be arbitrarily re-located as
$$F(t,ν) = \int f(t,p) δ(ν - ν(t,p)) ν d \log p$$
where the part of the scalogram at $(t,p)$ is relocated to $(t, ν(t,p))$, and still produces the same sound $f(t) = \int F(t,ν) d \log ν$. The actual relocation carried out is that corresponding to the Instantaneous Frequency, given here by the identity
$$4πi ν(t,p) |f(t,p)|^2 = f(t,p)^* \frac{∂}{∂t} f(t,p) - f(t,p) \frac{∂}{∂t} f(t,p)^*.$$
A simple, but naive, method (which was used in the video) is to just use the windowing $ψ(λ) = \exp(2πiλ)$ over one period $λ ∈ [-½,+½]$. It is sloppy, but gets the job done - mostly. However, with it, there will be effective interference, which registers as wiggling in the frequency lines in higher octaves. (That shows up in the video as the "wigglies" that are synched to frequencies below the displayed area - the sound is very deep. Use headphones!)
Windowing should actually be done in the frequency domain as
$$Ψ(γ) = |φ(\log γ)|^2, \quad ψ(λ) = \int |φ(\log γ)|² \exp(2πiγλ) dγ,$$
where $φ(z)$ is any of the usual spectrographic windowing functions, with a cut-off, say, of one octave $\log √½ ≤ z ≤ \log √2$, and normalized with $\int |φ(z)|^2 dz = 1$.
