I have a noisy signal x(t) that I want to differentiate in time, in order to obtain x'(t). Though numerical finite difference approximantion does not work well on noisy signal, I would like to differentiate its continuous wavelet transform (cwt) and to reconstruct the derivative from its wavelet transform. Is it possible to do so? Given the wavelet transform of the original signal, how do I differentiate it in order to get the wavelet transform of its derivative? I am using complex Morlet wavelet. Please bear in mind that I no matematician, so please be as human as you can in your answer.
[Math] Signal differentiation by wavelet transform
derivativeswavelets
Related Solutions
Well, since the wavelet transform is not precisely a time-frequency analysis but more a time-scale analysis one can not really display the wavelet coefficients in the time frequency plane. (However, even for the STFT the time-frequency plane is only a crutch for illustration purposes.)
Maybe it helps to think of $N=2^n$. Then you can analyze your signal into $n$ scales with the wavelet transform. See, e.g., this image with $n=3$:
Then the wavelet-transform stores the $c_{31}$ and all the $d_{ij}$'s. If you want to put the coefficients in the time-frequency plane you would have to stretch the tiles as you indicated in your figures. There would be $n$ vertical colums (in $\omega$), the rightmost having $2^{n/2}$ boxes and the number of boxes halving in each column (the last two containing just one box). This adds up to $2^n$.
However, I suggest to read Mallat's "A wavelet tour of signal processing - the sparse way" Chapters 1.3, 4.3 and 8.1. It's quite hard to produce illustrating figures for the wavelet transform here...
Considering the efficiency: To be a bit more precise, the discrete wavelet transform (DWT) of a signal of length $2^n$ and a filter with length $k$ takes $O(k2^n)$, while FFT is $O(n2^n)$ with comparable constants. Hence, DWT is faster it the filter is not too long.
Here's a proof of Morlet's inversion formula. We define \begin{align*} C_\psi &:= \int_0^\infty \dfrac{\overline{\hat{\psi}(u)}}{u}\, \mathrm{d}u \end{align*} and assume that this quantity is finite. Assume that $x,\psi \in L^2(\mathbb{R})$. I make this assumption so that Parseval's identity holds. Then, using Parseval's identity and basic properties of Fourier transforms (i.e., dilation becomes contraction on the Fourier side, and translation becomes phase modulation on the Fourier side), we can write
\begin{align*} X_w(a,b) &= \int_{-\infty}^\infty x(t) \frac{1}{\sqrt{a}} \overline{\psi\left(\frac{t-b}{a}\right)} \, \mathrm{d}t \\&= \int_{-\infty}^\infty \hat{x}(\xi) \frac{1}{\sqrt{a}} \overline{a\hat{\psi}\left(a\xi\right) \mathrm{e}^{-2\pi \mathrm{i} b \xi}}\, \mathrm{d}\xi \\&= \sqrt{a} \int_{-\infty}^\infty \hat{x}(\xi) \mathrm{e}^{2\pi \mathrm{i}b \xi}\overline{\hat{\psi}(a\xi)} \, \mathrm{d}\xi. \end{align*} Now, assume that $\hat{x} \in L^1(\mathbb{R})$ and that \begin{align*} A_\psi &:= \int_{0}^\infty \dfrac{\left|\hat{\psi}(u)\right|}{u} \, \mathrm{d}u \end{align*} is finite. I make these assumptions so that the precondition of the Fubini-Tonelli Theorem (i.e., absolute integrability of the double integral over $a$ and $\xi$) holds. Then, we can plug in our expression above for $X_w(a,t)$ and use the Fubini-Tonelli Theorem to switch the order of the integrals: \begin{align*} \frac{1}{C_\psi} \int_0^\infty X_w(a,t) \frac{\mathrm{d}a}{a^{3/2}} &= \frac{1}{C_\psi} \int_0^\infty \frac{\sqrt{a}}{a^{3/2}} \int_{-\infty}^\infty \hat{x}(\xi) \mathrm{e}^{2\pi \mathrm{i}b\xi} \overline{\hat{\psi}(a\xi)} \,\mathrm{d}\xi \mathrm{d}a \\&= \frac{1}{C_\psi} \int_{-\infty}^\infty \hat{x}(\xi) \mathrm{e}^{2\pi \mathrm{i}b\xi} \int_0^\infty \dfrac{\overline{\hat{\psi}(a\xi)}}{a} \, \mathrm{d}a \mathrm{d}\xi \\&= \frac{1}{C_\psi} \int_{-\infty}^\infty \hat{x}(\xi) \mathrm{e}^{2\pi \mathrm{i}b\xi} \int_0^\infty \dfrac{\overline{\hat{\psi}(u)}}{u} \, \mathrm{d}u \mathrm{d}\xi \\&= \frac{C_\psi}{C_\psi} \int_{-\infty}^\infty \hat{x}(\xi) \mathrm{e}^{2\pi \mathrm{i}b\xi} \\&= x(b), \end{align*} where in the third-to-last line we used the change-of-variable $u = \xi a$ in the inner integral, and in the last line we used Fourier inversion.
Note: The finiteness of $A_\psi$ is implied by the finiteness of $C_\psi$ provided that we're talking about Lebesgue integration. However, if we understand $C_\psi$ in the sense of an indefinite Riemann integral, then this is not necessarily guaranteed. (I can't come up with a real counterexample off the top of my head, but put for example $\hat{\psi}(\xi) = \sin(\xi)$; then $C_\psi$ is a finite indefinite Riemann integral but $A_\psi$ is infinite. Of course, this $\psi$ violates the assumption $\psi \in L^2$, so this isn't really a counterexample.)
Note 2: Regarding the idea of using different functions for reconstruction, the answer lies in the paper of Holschneider and Tchamitchian from 1991, Pointwise analysis of Riemann's "nondifferentiable" function (in Electricman's answer, the paper from Farge cites it as Holschneider & Tchamitchian 1989). Here, they prove the following theorem:
Let $\psi$ be as above and assume there exists a function $g \in L^1(\mathbb{R})$ where $g$ and $\psi$ together satisfy
- $a^{-1} \hat{g}(a) \hat{\psi}(a) \in L^1(\mathbb{R})$
- $\int_0^\infty a^{-1} \hat{g}(a) \overline{\hat{\psi}(a)}\, \mathrm{d}a = 1$
- $\int_0^\infty a^{-1} \hat{g}(-a) \overline{\hat{\psi}(-a)}\, \mathrm{d}a = 1$
- $\log(2 + |t|)\psi(t) \in L^1(\mathbb{R})$
- $\log(2 + |t|)g(t) \in L^1(\mathbb{R})$.
Suppose then that $x$ is a function which is "weakly oscillating about 0" in the sense that
$$ \lim_{u \rightarrow \infty} \dfrac{1}{2u} \int_{t-u}^{t+u} x(y) \, \mathrm{d}y = 0$$
where the limit holds uniformly in $t$.
Then at any point $t$ where $x$ is continuous, we have the inversion formula
$$ x(t) = \int_0^\infty \int_{-\infty}^\infty X_w(a,b) g\left(\dfrac{t-b}{a}\right)a^{-2} \, \mathrm{d}b \mathrm{d}a, $$
where the outer integral is understood to be a Riemann indefinite integral (i.e., you take limits as you approach the endpoints $0$ and $\infty$), and where $X_w(a,b)$ is the wavelet transform with respect to the mother wavelet $\psi$.
Of course, taking $g$ to be the delta function here isn't technically allowed (because $g$ must be in $L^1$) but it's still interesting to see this theorem.
Thank you for asking this really interesting question! I didn't know that such a different reconstruction could be done until I saw this post.
Best Answer
I usually work more with Discrete Wavelet Transforms, but I will give it a go. I base this guide below on an own rewritten version of Matlab source implementation of the Morlet wavelet.
$$ W(X) = ((\pi F_B)^{-0.5})\exp(2i\pi\cdot F_C \cdot X)\exp\left(-\frac{X^2}{F_B}\right)$$
It seems to me that the real and imaginary parts of the Morlet wavelet are in quadrature. This means they measure the whole phase range for a range of fourier frequencies. The first order derivative will be parts of the odd component, since a first order differential operator is, well odd. You have a bunch of parameters:
What is important if you want to create a differential approximator is to fine tune the bandwidth parameter. Too narrow band width and the wave will be very long, including many zero-crossings.
Example of Morlet wavelet where the wave is too long. We need to narrow it down in the time-domain, i.e. choose a shorter Gaussian window. The black envelope is the shape of our broad gaussian.
Here we can see that with a tighter Gaussian black envelope we pick out just one zero crossing for the odd part, which would be suitable to measure a derivative.