STFT and DWT (Wavelets)

Question

STFT can be successfully used on sound data (with a .wav soundfile for example) in order to do some frequency-domain modifications (example : noise removal).
With N=441000 (i.e. 10 seconds at sampling rate fs=44100), windowsize=4096, overlap=4, STFT produces approximatively a 430x4096 array (first coordinate : time frame, second coordinate : frequency bin). Modifications can be done on this array, and reconstruction can be done with overlap-add (*).

How is it possible to do a similar thing with wavelets ? (DWT), i.e. get a similar array of shape a x b, with a time frames, and b frequency bins, do some modification on this array, and at the end, recover a signal ? How ? What is the wavelet equivalent to overlap-add ? What would be the Python functions involved here (I haven't found an easy example of audio modification with pyWavelets...) ?

(*) : Here is the STFT framework that can be used :

signal = stft.Stft(x, 4096, 4)    # x is the input
modified_signal = np.zeros(signal.shape, dtype=np.complex)

for i in xrange(signal.shape[0]):    # Process each STFT frame
    modified_signal[i, :] =  signal[i, :] * .....  # here do something in order to
                                                   # modify the signal in frequency domain !
y = stft.OverlapAdd(modified_signal, 4)   # y is the output

The goal is to find a similar framework with wavelets.

Answer 1

The short-time Fourier transform is generally a redundant transformation, usually implemented with the same subsampling over every frequency. If the window is well chosen, it is complete: you can invert it and recover any initial signal.

Since it is redundant and complete, it has many perfect inverses. It can be implemented and understood using more generic tools: (oversampled) complex filter banks. Given a window type and length plus the overlap provides you with an analysis filter bank for which you can compute whether it is invertible or not. If it does, you can compute a natural inverse, and optimized inverses too. The overlap-add is just one on the many potential inverses, probably the most common, which often restricts the window choice.

Standard discrete wavelet transforms are filter banks too, with the difference that the subsampling is not the same in each frequency band (or more properly scale). This turns into uneven lengths for each scales. Yet, there exists redundant wavelet implementations that yield "a rectangle array" of coefficients you can work with. The best known schemes are called under different names: shift-invariant or time-invariant wavelets, undecimated wavelets, stationary wavelet transform (SWT), and sometimes cycle-spinning. Its standard reconstruction involves steps similar to overlap-add, except they are more "embedded" due to the different sampling factors over scales. You can use them with any discrete wavelet from a library, or even by designing your own wavelet. The reason is that standard discrete wavelets were designed with non-redundancy in mind, which restricts the choice of wavelets. With redundancy, you wavelet choice increases, as the constraints to fulfill are less stringent. The "ultimate" avatar is the continuous wavelet transform, which admits "almost" every inverse synthesis wavelet. My last sentence is quite lousy, I hope you get the meaning: when a square matrix is invertible, it has only one inverse. When a "rectangular" matrix is left-invertible in a generalized way, it has several inverses.

There seems to be a python implementation of the Stationary Wavelet Transform. You can find a few references in the 2.3.4. Translation invariant wavelets chapter of the linked paper.

It is generally much more robust for detection, denoising or restoration in practical applications (geophysics, non-destructive testing, ultrasounds, vibrations).

Answer 2

The reason you need overlap add / overlap save for filtering with the short time Fourier transform is basically, that the basis functions associated with the coefficients that you get are are defined over a certain time range (as opposed to a single point in time). The Fourier transform you use to calculate the expansion coefficients also implements convolution on a circular domain defined by your signal frame length. That means the two end points of the frame are really identified and closed in a circle. That is why you have to make sure that the basis functions of the coefficients you edit are never affecting both ends of the frame by wrapping around. This is a complication solely associated with discrete frame based Fourier transform methods and related to the convolution theorem and the fact that the discrete Fourier basis are the eigenvectors of time translations.

Wavelets are neither time translation eigenvectors nor are they calculated using circular convolution. This means you don't need overlap add or save or any other methods dealing with the side effects of circular convolution. Instead, the wavelet basis vectors are just a possible basis to describe your signal. The (complete, discrete, possibly orthogonal) wavelet transform is therefore nothing but a basis change from the time domain basis to the wavelet domain basis. Basis changes can be inverted (by applying the inverse of the basis change matrix that got you there) and you can switch back to time domain.

The parameters you gave as window size, overlap, sampling rate are all not applicable to the wavelet transform. The only thing you need is a mother wavelet. If you want to compare the results to your STFT output you can pick any of the STFT basis vectors (i.e your window multiplied by a complex exponential carrier) as the wavelet prototype. Then you apply the fast wavelet transform, which will decompose your signal into a tree of high and low pass filtered and decimated signals which ultimately become your coefficients. Each coefficient is associated with a wavelet basis vector and its parameters (scale,time) or (frequeny,time). You can manipulate the coefficients and then apply the inverse discrete wavelet transform. It will take your coefficients and run them through a resynthesis filter bank to produce a signal again.

These processes are not trivial and possibly hard to digest for a beginner. But you should be able to find libraries/toolboxes for your platform of choice that implement the fast wavelet transform and its inverse. However, if you want to realize your own wavelet basis, you will have to derive the filter coefficients for the decomposition and synthesis filter banks. That requires some deep theory, and you will probably have to study it first.

There are other flavors of the wavelet transform, namely the continuous wavelet transform which works with an overcomplete basis. It is both much slower to calculate and much harder to invert, so that it's not currently an option for what you want to do.

Answer 3

There are many ways of defining a wavelet basis. Typically a wavelet looks something like:

$$w_{x_0,k_0}(x) = A\exp(ik_0x)e(k_0(x-x_0))$$

In which $x_0$ is the centre in time, $k_0$ is the centre in frequency, and $e$ is a window function. $A$ absorbs the phase and normalisation. The main way in which this differs from your STFT is that the width of the window depends on $k$.

Usually one uses discrete points $(x_0,k_0)$ to limit the number of wavelets. As in your STFT, it is often good to use a larger number of points than the dimensionality of the signal. The aim is to approximate the answer you'd get by using all possible $(x_0,k_0)$, but with reasonable computational resources.

Because the dimensionality of the transformed data exceeds that of the signal, the wavelet basis will not be orthonormal. I.e. the following will be false:

$$\left\langle w_{k_0, x_0} | w_{k_0',x_0'} \right\rangle = \delta(x_0,x_0')\delta(k_0,k_0')$$

However, for suitable $A$ and $w$, you can arrange for the basis to be overcomplete:

$$\sum_{x_0,k_0} \left| w_{x_0,k_0} \right\rangle \frac{1}{\sqrt{k_0}} \left\langle w_{x_0,k_0} \right| = \text{identity}$$

In other words, you can reconstruct the signal perfectly just by adding up its constituent wavelets.

Your "modification" can simply be inserted in the above sum:

$$\text{my_filter} = \sum_{x_0,k_0} \left| w_{x_0,k_0} \right\rangle f(x_0,k_0) \left\langle w_{x_0,k_0} \right|$$

Update 2013-11-19: Adding implementation details below as requested.

For some signal $f(x)$, we wish to compute coefficients:

$$c_{x_0,k_0} = \left\langle w_{x_0,k_0} | f \right\rangle$$

For a fixed $k_0$, $c{x_0,k_0}$ may be viewed as a function of $x_0$, and this function is simply a filtered version of $f$. Specifically, it is a convolution of $f$ with $w_{0,k_0}$, which we can compute efficiently using a Fourier method. Thus, we may efficiently compute all the $c_{x_0,k_0}$ as follows:

• Apply a Fourier transform to $f$ to obtain $\hat f$. Probably you want to do this a window at a time, with enough overlap to be able to throw away the windowing artifacts, etc, but for simplicity let's suppose you do the whole signal at once, and that its length is a power of two.
• For each $k_0$ in some geometric progression with spacing about $1/4$ the filter bandwidth (or finer if you want):
  • Form the product of $\hat f$ with $\hat w_{0,k_0}$.
  • Truncate the spectrum to some interval $[k_l,k_r)$ whose length is a power of two, and that contains the non-zero portion of $\hat w_{0,k_0}$.
  • Apply an inverse Fourier transform to that.
  • Multiply that by $exp(ix\frac{k_l+k_r}{2})$ to correct the phase. The result is $c_{x_0,k_0}$ viewed as a function of $x_0$.

This computes all the wavelet coefficients. You can choose the resolution in $k_0$ by tuning the ratio in the geometric progression. The resolution in $x_0$ is set by the length of the truncated spectrum, and will change depending on the bandwidth of $w_{0,k_0}$, which in turn depends on $k_0$. The computational effort is one Fourier transform at high time resolution, plus one inverse Fourier transform for each $k_0$ value at much lower time resolution. It works out about the same as STFT - maybe slower by some small factor that depends on the resolution you choose.

You can then modify the $c_{x_0,k_0}$ as you see fit, and you can reconstruct the signal by reversing the above process, summing the spectra over $k_0$ before finally doing an overall inverse Fourier transform.

Truncating spectra sometimes introduces normalisation problems, depending on precisely how your FFT is defined. I will not attempt to cover all the possibilities here. Normalisation is basically an easy problem. ;-)

The only part that remains is to choose a suitable wavelet envelope. It turns out it is easier to get $\hat w_{x_0,k_0}(k)$ right than to get $w_{x_0,k_0}(x)$ right. A suitable definition (from many possibilities) is:

$$\hat w_{x0,k0} = A exp(-i(k-k_0)x_0) exp(-(Q log(k/k_0))^2)$$

in which $Q$ is a dimensionless constant that selects the bandwidth of your filter, i.e. the frequency resolution of your wavelets, and $A$ is chosen as necessary for normalisation. With this definition and sufficiently high resolution for $k_0$, the overcompleteness condition holds, and the signal reconstruction will work.