Advantage of STFT over wavelet transform

Question

I have learned about STFT and wavelet transform recently, and wavelet transform seems better than STFT in my opinion. So, I wonder if there is any advantage of using STFT than WT, and if so, what are practical applications of STFT?

Answer 1

Wavelet transforms and short-term/short-time Fourier transforms are broad names for classes of transformations that are not totally distinct and may overlap (pun intended).

Both can be efficient for non-stationary features of data, and they both have merits or drawbacks, depending on their parameters and signal's properties. STFT is typically analyzing signals on fixed-length windows with different modulations, while wavelets are similar modulations (zero-crossing) on different support sizes.

I am a promoter of wavelet-type methods. I however should mention that in image and audio JPEG and mp3 are widely spread standards akin the STFT (fixed length), in their critical version (maximally sub-sampled). Wavelets, although at the base of the JPEG2000, is less used, possibly for implementation/use issues.

In video coding, and in deep learning, it is more customary to look at different resolutions (akin to wavelets), yet not exactly in the structured dyadic wavelet fashion.

Answer 2

STFT is frequency-shift equivariant - same absolute shift has same effect on representation regardless of original frequency${}^1$:

$$ \hat x(\omega) \rightarrow \hat x(\omega - c) \Leftrightarrow \text{STFT}_x (t, \omega) \rightarrow \text{STFT}_x (t, \omega - c) \tag{1} $$ _{1: rather $|\omega| - c$ in general, if considering negative freqs}

This allows it to track linear frequency modulation with a fixed rate over the entire time-frequency plane, so we can tune the time-frequency tradeoff until we've mapped it perfectly. Shift equivariance is also useful for representative linearity depending on application (e.g. quantifying audio frequency transposition).

An overlooked advantage is that, STFT is much easier to implement - even major Python libraries (PyWavelets, scipy) have flaws. It can also be considered faster, per more permissive "hop size".

Overall I do favor CWT over STFT - with CWT properties in depth here (rather scattering, but some apply to CWT also). Lastly, more comparisons on test signals.

Answer 3

In response to the reply on unifying the STFT and CWT (I can't comment yet): I keep a recoded version of ARSS over here

https://github.com/LydiaMarieWilliamson/ARSS

though it will undergo rebasing and some revisions soon - partly to incorporate the new ideas. The scaling isn't actually done quite in that way in ARSS, but a little differently. Instead of going like $P = \left(e^{kx} - 1\right)/k$, it's better to do it like this: $P = \left(e^{kx} - 1\right)/\left(e^k - 1\right)$, then you can say that the case $k = ln(2)$ produces a simple exponential, $P = 2^x - 1$. This shows that there should also be an offset - which is what ARSS actually has in it, so that the general form would be $P(x) = P₀ B^x + β \left(B^x - 1\right)/\left(B - 1\right)$ - which is how it's actually done by ARSS (up to a change in notation and names). With this scaling, we would then have $ΔP(x) = P(x + 1) - P(x) = α P(x) + β$ with $α = B - 1$. So, the respective cases would be $α = 0$ for the time frequency plane - corresponding to the limit $B → 1$, where $P(x) → P₀ + βx$ and $β = 0$ for the time-scale plane, with the limit $P(x) = P₀ B^x$.

So, that's the generalization that is linear for small frequencies and exponential for high frequencies.

Missing in your description is what Lie algebraic underpinnings there would be to the unified formulation. Since the time-frequency plane corresponds to the Lie algebra of the Heisenberg group, and the time-scale plane to the dilation group (with respective Lie brackets $[q,p] = z$, $[q,z] = 0$, $[p,z] = 0$ in the former case, and $[d,p] = p$ in the latter case), then what corresponds to the unified generalization of the two? If we take your suggestion literally, then perhaps we should first throw together all the generators (and add the Lie brackets $[q,d] = q$ and $[d,z] = 0$) and then restrict to the subalgebra with $y = αd + βq$, $p$ and $z$, to yield the Lie algebra $[y,p] = αp + βz$, $[y,z] = 0$ and $[p,z] = 0$. Then the case $\left(α,β\right) = \left(1,0\right)$ yields the Lie algebra of the dilation group (along with a trivial central extension of it by $z$), while the case $\left(α,β\right) = \left(0,1\right)$ yields the Lie algebra of the Heisenberg group.

Also left unanswered is: how could the CWT, or (more generally) its unified generalization, be streamed? Since an infinite wave train has infinite bandwidth, then the base frequency is infinitesimal. So, the transform would either have to go linear for small frequencies, so that you could get down to zero, in which case $β ≠ 0$ is necessary.

Otherwise, the stream would have to be chopped up - which (however) may put restrictions on what kinds of windowing functions $h$ you could use for the inverse, if you also want a single unified graph for the entire stream that admits an inverse that yields the infinite stream. You might need something like an STFT stage, with really large step sizes (e.g. 1 second), and then run the transform on each (1-second) window. This might be a related reference in that vein

Xisheng Li and Shaochun Wang, "A WT-STFT combining Algorithm", J. Univ. Sci. Technol. Beijing, 7(2000), No. 4, pp. 315-317. http://ijmmm.ustb.edu.cn/en/article/id/947345fd-718c-4d41-a471-87d01cdf9e16

It doesn't appear to actually be talking about hybridizing the two tranform families, as you just did, but actually to be using an STFT stage for a CWT transform, along the lines of what I just suggested. Perhaps, this could be generalized to the entire "CWFT" family, as you call it.

Answer 4

First, a quick answer, then a lengthy reply calling into question the question, itself. The main advantage of the windowed Fourier transform (or STFT, as you call it) is that it can be done cleanly on an infinite stream with finite resources - especially because of the regularity of its underlying grid. You can do it in real time on streaming signals.

Though the time increment $Δt$ arises naturally in the signal (as the sampling period), the corresponding frequency increment $Δν$ does not arise from the signal, as it would be $Δν = 1/T$, for a signal of finite duration $T$, if it did, and with the signal being unlimited in length $T → ∞$, then $Δν → dν$ an infinitesimal and the corresponding Fourier transform would be 1/2 continuous and 1/2 discrete. However, a finite $Δν$ does arise from the width of the windowing function used with the transform. The window cuts off the signal into finite-duration chunks, each one which is converted in the usual way.

In contrast, time-scale transforms make essential use of the frequency domain; and it would be very difficult to make it work with streaming signals; particularly, as the sampling points in the time-scale plane are irregularly-spaced - which is normally touted as a strength. Here, though, the strength of the transform is a weakness when it comes to the issue of discretization.

It's been an active area of research; example:

Peter Heinlein, "Discretizing continuous wavelet transforms using integrated wavelets"
Applied and Computational Harmonic Analysis, Volume 14, Issue 3, May 2003, Pages 238-256
https://www.sciencedirect.com/science/article/pii/S1063520303000058

A time-scale transform could be done in the time-domain, though it's not clean. One transform which I devised (actually: stumbled onto, by accident) is $$\tilde f(q,p) = \int_{-1/2}^{+1/2} f\left(q + {λ \over p}\right) e^{-2πiλ} dλ.$$ The inverse is just simple summation on a logarithmic scale, since: $$∫ \tilde f(q,p) δ(t-q) {dq dp \over |p|} = ∫ \tilde f(t,p) {dp \over |p|} = f(t),$$ so the time series for $q$, of each $p$ channel (corresponding to the scale $1/p$) is a bona-fide signal component in its own right. The wavelet, here, for the forward transform is just $ψ_{q,p}(t) = exp(2πip(t-q)) |p|$ cut off to a single period in the range $t ∈ [q-1/2p,q+1/2p]$, while for the inverse transform it's not even a windowing function, at all, but just a spike: $δ(t-q)/|p|$.

Either-Or?
Instead of asking this question, you should first be asking the question: is it even an either-or situation, in the first place? That is: is it possible to actually unify the wavelet and windowed Fourier transform seamlessly into a single framework, in such a way that each becomes special cases of the more general transform?

I first encountered this issue when reviewing the code for the open source ARSS program (Analysis and Resynthesis Sound Spectrograph).

http://arss.sourceforge.net/

It's abandonware now, as the author has gone on to make a commercial replacement for it that has a GUI front-end and a broader range of functions.

The author, at the time, was not aware that he actually devised a slightly generalized form of the wavelet transform. But one of the most interesting features of the software is that it has a parameter which, when set to 0, yielded a linear scale suitable for the windowed Fourier transform, and when set to a non-zero cardinal value (the natural log of 2) yielded a logarithmic scale suitable for the wavelet transform.

Unfortunately, he was never able to generalize it beyond the two cardinal values, so the software effectively became a disjoint union of the two transforms, rather than a bridging between the two. The same is true of the commercial version.

So, the question was left unresolved: can the gap between the two actually be bridged? And the answer is - if you make the right generalizations to each of the two families of transforms (with the generalization required for the wavelet transforms also resulting in bringing the S transforms in as special cases), you can actually devise a unified framework that seamlessly combines the two. I'll describe that here, in more detail.

As they are normally defined and formulated, they cannot be combined, because they are defined in a way that imposes unnecessary conditions - and the respective definitions are more constrained than they need to be. So, first, let's look at the respective transforms and make the appropriate repairs to each of them, before going on.

For the following, to simplify the notation (tremendously!), I will be adopting the following convention: $$1^x ≡ e^{2πix} = \cos 2πx + i \sin 2πx.$$ Under this convention, the Fourier transform of a function $f(t)$ for $t ∈ ℝ$, and its inverse are given, respectively, by: $$\hat f(ν) = ∫ f(t) 1^{-νt} dt,$$ $$f(t) = ∫ \hat f(ν) 1^{νt} dν,$$ where $ν$ is the linear frequency - which (as you'll see) makes things much simpler than if using the angular frequency $ω = ν/{2π}$. Combining the two transforms yields the respective identity $$f(t) = ∫∫ f(t') 1^{ν(t'-t)} dt' dν = ∫ f(t') \left(∫ 1^{ν(t'-t)} dν\right) dt',$$ thus showing that the kernel of the transform is just the delta function $$δ(T) = ∫ 1^{νT} dν,$$ which is defined by the condition that $$f(t) = ∫ f(t') δ(t'-t) dt'.$$

Windowed Fourier Transform (WFT)
Also called the Short-Time Fourier Transform (STFT), the WFT is a windowed transform that uses a windowing function of the following form: $$ψ_{q,p}(t) = g(t-q) 1^{p(t-q)},$$ where $q$ and $p$ are, respectively the horizontal and vertical coordinates, suitable for the representation on a spectrogram of a time series given by the sequences of $q$'s for multiple frequency channels given by the different values of the $p$'s.

This is a complex sinusoidal modulated by an envelope, given by the function $g$. The usual shape for $g(Q)$ is a hump centered around $Q = 0$ with rapidly dropping tails as you get away from $0$. Nothing, either above or in the following, requires $g(Q)$ to actually be real-valued, though that's the case for all the choices of envelope functions $g$ that I know of.

Under the transform, a function $f(t)$ of one coordinate (the time $t$) is converted to a function over the Time-Frequency Plane (with coordinates $q,p$), given by $$\tilde f(q,p) ≡ ∫ f(t) \overline{ψ_{q,p}(t)} dt$$ where $\overline{\left(⋯\right)}$ denotes complex conjugation.

Under certain, very broad, conditions the transform can be inverted - with the inverse given as a windowed transform over the time-frequency plane, itself. To formulate those conditions, suppose we apply the transform with the following choice for the windowing function $$φ_{q,p}(t) = h(t-q) 1^{p(t-q)}.$$ Then, the condition we want is for the following integral $$F(t) ≡ ∫ \tilde f(q,p) φ_{q,p}(t) dq dp$$ to yield $f(t)$, itself.

Normally, for the WFT (as best as I understand it), the same windowing function is used in both directions, so that $h = g$ and $φ = ψ$. But $h = 1$ could just as well be used. It's actually more convenient, because with that choice, the spectrogram is showing you actual wave forms for the components of the time series given by the function $f(t)$.

It is also common to insert another factor of $1^{pq}$ is the forward transform, so that the windowing function would actually be $ψ_{q,p}(t) = g(t-q) 1^{pt}$, instead; and a $1^{-pq}$ in the inverse transform to make the windowing function $φ_{q,p}(t) = h(t-q) 1^{pt}$. The values on the spectrogram are complex. One way to denote them is by displaying the amplitude by brightness and the phase by color. However, beyond a certain resolution you won't see the color (it will all be whited out) in a graph. Without the insertion $1^{pq}$, the phase, at each level $p$ is the actual phase of the waveform at $p$ - particularly, if the function $h(Q) = 1$ is chosen for the inverse transforms. With the inserted factor $1^{pq}$, the phase tends toward a constant at those values of $p$ that best match the actual underlying frequency of the component.

To find what the conditions on $g$ and $h$ are, it is best to write the transform and its inverse in terms of the frequency coordinate $ν$, instead of $t$. This is done, here, first in general form in terms of the windowing functions $ψ_{q,p}$ and $φ_{q,p}$, assuming first only the following form of translational invariance: $$ψ_{q,p}(t) = ψ_p(t - q),$$ $$φ_{q,p}(t) = φ_p(t - q).$$ For the forward transform, we have $$\tilde f(q,p) = ∫∫∫ \hat f(ν) 1^{νt} \overline{\left(\hat {ψ_p}(ν') 1^{ν'(t-q)}\right)} dν dν' dt$$ The $t$'s are all up in the exponentials, so can be factored and integrated out to a delta function $$∫ 1^{νt} \overline{1^{ν't}} dt = ∫ 1^{(ν-ν')t} dt = δ(ν-ν').$$ So, upon substitution, and integration over $ν'$ (which becomes trivial, because of the delta function) we get $$\tilde f(q,p) = ∫∫ \hat f(ν) δ(ν-ν') \overline{\hat {ψ_p}(ν')} 1^{ν'q} dν dν' = ∫ \hat f(ν) \overline{\hat {ψ_p}(ν)} 1^{νq} dν.$$

For the inverse direction, in the frequency domain, we have $$\hat F(ν) = ∫∫∫∫ \tilde f(q,p) 1^{-νt} \hat{φ_p}(ν') 1^{ν'(t-q)} dt dν' dq dp.$$ Again, the $t$'s are all up in the exponential, so by a similar process the factor $1^{(ν'-ν)t}$ can be collected, integrated out by $t$ to a delta function $δ(ν'-ν)$, and then the $ν'$ integrated out (again, becoming trivial because of the delta function) to yield: $$\hat F(ν) = ∫∫ \tilde f(q,p) \hat{φ_p}(ν) 1^{-νq} dq dp.$$ When combined with the expression for $\tilde f(q,p)$, this results in: $$\hat F(ν) = ∫∫∫ \left(\hat f(ν') \overline{\hat {ψ_p}(ν')} 1^{ν'q}\right) \hat{φ_p}(ν) 1^{-νq} dν'dq dp,$$ and now (because of the assumption made about translational invariance for $ψ_{q,p}$ and $φ_{q,p}$), the $q$'s are all up in the exponentials, the factors can be collected to $1^{(ν'-ν)q}$, integrated with respect to $q$ to yield the delta function $δ(ν'-ν)$, as before, with $ν'$ integrated out to yield: $$\hat F(ν) = \hat f(ν) ∫ \overline{\hat {ψ_p}(ν)} \hat{φ_p}(ν) dp = \hat f(ν) \hat Γ(ν).$$ So, without further conditions on the windowing functions, the result is the same as doing a filter with the function $$\hat Γ(ν) = ∫ \overline{\hat {ψ_p}(ν)} \hat{φ_p}(ν) dp,$$ with $F(t)$ being the convolution of $f(t)$ and the time-domain version $Γ(t)$ of $\hat Γ(ν)$.

When expressed in terms of the envelope functions $g$ and $h$, we have $$\hat {ψ_p}(ν) = ∫ g(Q) 1^{pQ} 1^{-νQ} dQ = \hat g(ν-p),$$ $$\hat {φ_p}(ν) = ∫ h(Q) 1^{pQ} 1^{-νQ} dQ = \hat h(ν-p),$$ and, finally, we find that $$\hat Γ(ν) = ∫ \overline{\hat g(ν-p)} \hat h(ν-p) dp = ∫ \overline{\hat g(P)} \hat h(P) dP = \left<\hat g,\hat h\right> = \left<g,h\right> = ∫ \overline{g(Q)} h(Q) dQ.$$ For the choice $h(Q) = 1$ (i.e. $\hat h(P) = δ(P)$), this reduces to the condition $\hat g(0) = 1$. For $h$ = $g$, it is just the condition that ${|g|}^2 = \left<g,g\right> = 1$. In general, the requirement that $F = f$ is that $\left<g,h\right> = 1$.

Continuous Wavelet Transform (CWT)
The wavelet transform is a windowed transform which uses a windowing function of the form $$ψ_{στ}(t) = ψ\left({t-τ}\over σ\right) {1\over \sqrt{|σ|}},$$ where $τ$ and $σ$ are, respectively, the time and scale in the Time-Scale Plane. We can just as well write this as just another case of a windowing function already treated in general, by the correspondence $$σ ⇔ {1 \over p}, τ ⇔ q$$ and write $$ψ_{q,p}(t) = ψ(p(t-q)) {|p|}^A,$$ with $A = ½$.

(I will ignore the issue, in the following, of what to do about 0 or negative values of $p$.)

The inverse transform is normally taken with the same windowing function as the forward transform. The restrictions on $ψ$ which result from this are severe enough to prevent natural windowing functions from being used - like Gaussian-enveloped complex sinusoidal functions. One can only obtain $\hat Γ(ν) = 1$ for certain, very kooky and arcane shaped functions $ψ$.

However, there is nothing preventing us from generalizing by allowing different windowing functions to be used in the forward and reverse directions. So, we will write the windowing function for the reverse direction as $$φ_{q,p}(t) = φ(p(t-q)) {|p|}^B.$$ (Note: this notation conflicts with the convention, in the wavelet literature, of using $φ$ to denote another function, which is unrelated to the one we defined here.) In the case where the two sets of windowing functions are the same, we have $φ = ψ$ and $B = ½$.

Again, as before, we have a form of translational invariance: $$ψ_{q,p}(t) = ψ_p(t-q), φ_{q,p}(t) = φ_p(t-q),$$ where $$ψ_p(Q) = ψ(pQ) {|p|}^A, φ_p(Q) = φ(pQ) |p|^B.$$ The respective Fourier transforms are $$\hat{ψ_p}(ν) = ∫ ψ(pQ) 1^{-νQ} {|p|}^A dQ = ∫ ψ(T) 1^{-{νT \over p}} {|p|}^{A-1} dT = \hat ψ\left({ν \over p}\right) {|p|}^{A-1},$$ and similarly, $$\hat{φ_p}(ν) = ∫ φ(pQ) 1^{-νQ} {|p|}^B dQ = ∫ φ(T) 1^{-{νT \over p}} {|p|}^{B-1} dT = \hat φ\left({ν \over p}\right) {|p|}^{B-1}.$$

Upon substitution, the condition we find that $$\hat Γ(ν) = ∫ \overline{{\hat ψ}\left({ν \over p}\right)} {\hat φ}\left({ν \over p}\right) {|p|}^{A+B-1} {dp \over |p|}.$$ Assuming that the domains of ψ and φ are restricted only to positive values, then after making the substitution $ν/p = e^x$, this reduces to $$\hat Γ(ν) = ∫ \overline{{\hat ψ}\left(e^x\right)} {\hat φ}\left(e^x\right) {\left|ν e^{-x}\right|}^{A+B-1} dx.$$

The requirement that $\hat Γ(ν)$ be independent of $ν$ means that we must have $A + B = 1$, with the standard choice being $A = ½ = B$, but actually $A = 0$ and $B = 1$ being more convenient. Under this assumption, the condition that $\hat Γ(ν) = 1$ yields the result that $$∫ \overline{{\hat ψ}\left(e^x\right)} {\hat φ}\left(e^x\right) dx = 1,$$ which shows that the frequencies are naturally on a logarithmic scale, in the time-scale plane, with $x = \ln(ν/p)$ being equally-spaced between octaves.

This can be just as well written as enveloped sinusoidals by defining $g(Q) = ψ(Q) 1^{-Q}$ and $h(Q) = φ(Q) 1^{-Q}$ (i.e. $\hat g(P) = \hat ψ(P+1)$ and $\hat h(P) = \hat φ(P+1)$). Then, the corresponding condition on $g$ and $h$ become $$1 = ∫ \overline{{\hat g}\left(P\right)} {\hat h}\left(P\right) {dP \over |P+1|} = ∫ \overline{{\hat g}\left(e^x - 1\right)} {\hat h}\left(e^x - 1\right) dx.$$ There is lots of room available to permit $ψ$ to be just a sinusoidal enveloped by $g$, like with the WFT, e.g. just taking $\hat g(P)$ to be a cosine window in $x$, itself, with $\hat h(P) = 1$. This has close connection to the so-called S-transform (it's a generalization of it) and, I think, also to the Q transform (which might correspond to taking $h(Q) = 1$, but I'm not sure).

CWFT: Unification of the WFT and CWT
A short while back, I posted a reply - elsewhere - that Gabor ... who had devised the WFT in 1946 ... missed the boat on the CWT, because he actually had a little-known section in his pioneering paper discussing the nature of human sound perception, which noted the logarithmic scale, but made an unwarranted assumption in his paper that the windowing functions should correspond to constant window sizes $Δp$ and $Δq$. For the CWT, the window sizes go as $Δp$ proportional to $p$, while $Δq$ is proportional to $1/p$ or to the scale $σ$.

Actually, that was partly a mistaken attribution. I think it was Daubechies who pointed out the details of the mechanisms for hearing pitch on a logarithmic scale in her pioneering work on the wavelet transform, rather than Gabor. The paper by Gabor does discuss human sound perception, but actually points out that the natural windowing size for $Δp$ is constant up to around 1000 Hz, and only then starts to become logarithmic for higher frequencies, with $Δp$ becoming proportional to $p$. In other words, the scaling may actually be hybrid.

This may be what the author of ARSS was trying to embody.

Now that they've been brought into a common form, it is actually possible to take it one step further and write this as two cardinal points of an entire family of transforms, with the two extremes bridged by the cases of hybrid transforms.

First, in order for them to be unified, we have to treat the enveloping functions $g$ and $h$, independently of the sinusoidal in the windowing functions $ψ$ and $φ$, since they were using different combinations of $p$, $q$ and $t$ in the WFT than what appeared in the complex sinusoidals. In the CWT, the same combination $p(t-q)$ appears in both the enveloping functions and in the complex sinusoidal - and that's what needs to be generalized away from.

The ARSS software use a hybrid linear-exponential scale of the form $a + b e^{kx}$ that reduced a linear progression in $x$ as $k → 0$. With that in mind, we first generalize the integral expression for $\hat Γ(ν)$ in terms of $x$ by replacing the $x$ of the WFT and the $e^x - 1$ of the CWF by $\left(e^{kx} - 1\right)/k$: $$\hat Γ(ν) = ∫ \overline{\hat g\left({e^{kx} - 1}\over k\right)} \hat h\left({e^{kx} - 1}\over k\right) dx,$$ or, upon substitution of $P = \left({e^{kx} - 1}\right)/k$, $$\hat Γ(ν) = ∫ \overline{\hat g\left(P\right)} \hat h\left(P\right) {dP \over \left|kP+1\right|}.$$

Second, we can then combine the forms $P = ν - p$ for the WFT and $P = \left(ν - p\right)/p$ for the CWT into $P = \left(ν - p\right)/\left(kp + l\right)$, with the WFT and CWT being the respective special cases of $(k,l) = (0,1)$ and $(k,l) = (1,0)$. Then, we may write: $$\hat Γ(ν) = ∫ \overline{\hat g\left({ν - p \over kp + l}\right)} \hat h\left({ν - p \over kp + l}\right) {dp \over \left|kp + l\right|}.$$

From this, we can read off the respective transforms of the windowing functions which - upon restoring the integrating factors with $A + B = 1$, as before - are: $$\hat{ψ_p}(ν) = g\left({ν - p}\over{kp + l}\right) \left|kp + l\right|^{A-1},$$ $$\hat{φ_p}(ν) = h\left({ν - p}\over{kp + l}\right) \left|kp + l\right|^{B-1},$$ or in back in the time domain as: $$ψ_p(Q) = g\left((kp + l)Q\right) 1^{pQ} \left|kp + l\right|^A,$$ $$φ_p(Q) = h\left((kp + l)Q\right) 1^{pQ} \left|kp + l\right|^B,$$ or, finally, substituting back $Q = t - q$: $$ψ_{q,p}(t) = g\left((kp + l)(t - q)\right) 1^{p(t - q)} \left|kp + l\right|^A,$$ $$φ_{q,p}(t) = h\left((kp + l)(t - q)\right) 1^{p(t - q)} \left|kp + l\right|^B.$$

Related References
I'm not the only one to consider these issues. A recent reference:

Carlos Mateo & Juan Antonio Talavera
Bridging the gap between the short-time Fourier transform (STFT), wavelets, the constant-Q transform and multi-resolution STFT
Signal, Image and Video Processing, volume 14, pages 1535–1543 (2020)
https://link.springer.com/article/10.1007/s11760-020-01701-8

and another from further back:

Xisheng Li and Shaochun Wang, A WT-STFT combining Algorithm,
J. Univ. Sci. Technol. Beijing,
International Journal of Minerals, Metallurgy and Materials 2000 7(4) 315-317
http://ijmmm.ustb.edu.cn/article/id/947345fd-718c-4d41-a471-87d01cdf9e16