# Daubechies wavelet transform

i have N samples obtained by sampling a signal with lot of frequency contents. How will i apply daubechies wavelet transform to obtain the frequency and its location? i need to write a program which will process the signal and gives the frequency and location as the result.

• If you have multiple frequencies (more like bursts of frequency?) in the signal, which one are you interested of? – Olli Niemitalo Feb 3 '16 at 7:20
• i am interested in getting the higher frequencies(i.e above a threshold which i set) present. – vvv Feb 3 '16 at 8:34
• Higher as in higher frequency or higher amplitude? – Olli Niemitalo Feb 3 '16 at 8:53
• if possible i would like to know the amplitude of the signal at those higher frequencies as well. – vvv Feb 3 '16 at 8:59
• Have you managed to do the transform already? – Olli Niemitalo Feb 3 '16 at 9:13

Looks like you need a general explanation of the discrete wavelet transform (DWT). DWT breaks a signal down into subbands distributed evenly in a logarithmic frequency scale, each subband sampled at a rate proportional to the frequencies in that band. The traditional Fourier transformation has no time domain resolution at all, or when done using many short windows on a longer data, equal resolution at all frequencies. The distribution of samples in the time and frequency domain by DWT is of form:

log f
|XXXXXXXXXXXXXXXX  X = a sample
|X X X X X X X X   f = frequency
|X   X   X   X     t = time
|X       X
|X
|X
----------------t


Single subband decomposition and reconstruction:

    -> high -> decimate -------------> dilute -> high
|   pass    by 2      high subband  by 2      pass \
in |                                                   + out
|                                                  /  =in
-> low  -> decimate -------------> dilute -> low
pass    by 2      low subband   by 2      pass


This creates two subbands from the input signal, both sampled at half the original frequency. The filters approximate halfband finite impulse response (FIR) filters and are determined by the choice of wavelet. Using Daubechies wavelets (and most others), the data can be reconstructed to the exact original even when the halfband filters are not perfect. Note that in the above scheme, the total amount of information (samples) stays the same throughout.

Decimation by 2: ABCDEFGHIJKLMNOPQR -> ACEGIKMOQ
Dilution by 2: ACEGIKMOQ -> A0C0E0G0I0K0M0O0Q0


To get the logarithmic resolution in frequency, the low subband is re-transformed, and again, the low subband from this transformation gets the same treatment etc.

Decomposition:

    -> high -> decimate --------------------------------> subband0
|   pass    by 2
in |                 -> high -> decimate ---------------> subband1
|                |   pass    by 2
-> low -> decim |                 -> high -> decim -> subband2
pass   by 2  |                |   pass    by 2
-> low -> decim |
pass   by 2  |   .   down to what suffices
-> .    or if periodic data,
.     until short of data


Reconstruction:

subband0 -----------------------------------> dilute -> high
by 2      pass \
subband1 ------------------> dilute -> high                   + out
by 2      pass \                /  =in
subband2 -> dilute -> high                   + dilute -> low
by 2      pass \                /  by 2      pass
+ dilute -> low
Start   .                  /  by 2      pass
here!   . -> dilute -> low
.    by 2      pass


In a real-time application, the filters introduce delays, so you need to compensate them by adding additional delays to less-delayed higher bands, to get the summation work as intended.

For periodic signals or windowed operation, this problem doesn't exist - a single subband transformation is a matrix multiplication, with wrapping implemented in the matrix:

Decomposition:

|L0|   |C0  C1  C2  C3                | |I0|     L = lowpass output
|H0|   |C3 -C2  C1 -C0                | |I1|     H = highpass output
|L1|   |        C0  C1  C2  C3        | |I2|     I = input
|H1| = |        C3 -C2  C1 -C0        | |I3|     C = coefficients
|L2|   |                C0  C1  C2  C3| |I4|
|H2|   |                C3 -C2  C1 -C0| |I5|
|L3|   |C2  C3                  C0  C1| |I6|
|H3|   |C1 -C0                  C3 -C2| |I7|     Daubechies 4-coef:

1+sqrt(3)        3+sqrt(3)        3-sqrt(3)        1-sqrt(3)
C0 = ---------   C1 = ---------   C2 = ---------   C3 = ---------
4 sqrt(2)        4 sqrt(2)        4 sqrt(2)        4 sqrt(2)


Reconstruction:

|I0|   |C0  C3                  C2  C1| |L0|
|I1|   |C1 -C2                  C3 -C0| |H0|
|I2|   |C2  C1  C0  C3                | |L1|
|I3| = |C3 -C0  C1 -C2                | |H1|
|I4|   |        C2  C1  C0  C3        | |L2|
|I5|   |        C3 -C0  C1 -C2        | |H2|
|I6|   |                C2  C1  C0  C3| |L3|
|I7|   |                C3 -C0  C1 -C2| |H3|


C0, C1, C2, C3 are the "db2" lowpass FIR filter coefficients. Highpass coefficients you get by reversing tap order and multiplying by sequence 1,-1, 1,-1, ... Because these are orthogonal wavelets, the analysis and reconstruction coefficients are the same.

A coefficient set convolved by its reverse is an ideal halfband lowpass filter multiplied by a symmetric windowing function. This creates the kind of symmetry in the frequency domain that enables aliasing-free reconstruction. Daubechies wavelets are the minimum-phase, minimum number of taps solutions for a number of vanishing moments (seven in "db7" etc.), which determines their frequency selectivity.

I was asked to show the matrices for 6 coefficients, so here they are, made a bit larger for clarity but could be the same size as before too. Decomposition:

|L0|   |C0  C1  C2  C3  C4  C5                | |I0|
|H0|   |C5 -C4  C3 -C2  C1 -C0                | |I1|
|L1|   |        C0  C1  C2  C3  C4  C5        | |I2|
|H1|   |        C5 -C4  C3 -C2  C1 -C0        | |I3|
|L2| = |                C0  C1  C2  C3  C4  C5| |I4|
|H2|   |                C5 -C4  C3 -C2  C1 -C0| |I5|
|L3|   |C4  C5                  C0  C1  C2  C3| |I6|
|H3|   |C1 -C0                  C5 -C4  C3 -C2| |I7|
|L4|   |C2  C3  C4  C5                  C0  C1| |I8|
|H4|   |C3 -C2  C1 -C0                  C5 -C4| |I9|


Reconstruction:

|I0|   |C0  C5                  C4  C1  C2  C3| |L0|
|I1|   |C1 -C4                  C5 -C0  C3 -C2| |H0|
|I2|   |C2  C3  C0  C5                  C4  C1| |L1|
|I3|   |C3 -C2  C1 -C4                  C5 -C0| |H1|
|I4| = |C4  C1  C2  C3  C0  C5                | |L2|
|I5|   |C5 -C0  C3 -C2  C1 -C4                | |H2|
|I6|   |        C4  C1  C2  C3  C0  C5        | |L3|
|I7|   |        C5 -C0  C3 -C2  C1 -C4        | |H3|
|I8|   |                C4  C1  C2  C3  C0  C5| |L4|
|I9|   |                C5 -C0  C3 -C2  C1 -C4| |H4|


With:

C0 = 3.326705529500826159985115891390056300129233992450683597084705e-01
C1 = 8.068915093110925764944936040887134905192973949948236181650920e-01
C2 = 4.598775021184915700951519421476167208081101774314923066433867e-01
C3 = -1.350110200102545886963899066993744805622198452237811919756862e-01
C4 = -8.544127388202666169281916918177331153619763898808662976351748e-02
C5 = 3.522629188570953660274066471551002932775838791743161039893406e-02


More coefficient sets can be found here.

• Thanks a lot sir, i understood about this decomposition and reconstruction, but i am unable to understand ,like if i write a code on this dwt, how will i decipher the frequency and time details. More clearly, i wanted to do like, i am giving 32 samples as input and my code should give the frequency,time content of the samples whose amplitude exceeds a specific threshold. – vvv Feb 3 '16 at 12:01
• In the DWT results, the subband index of a sample tells the frequency, and its time index together with the subband index (which determines the sampling frequency of the subband) tells the time. – Olli Niemitalo Feb 3 '16 at 15:02
• Sir, i have understood like, to decompose and reconstruct the signals, we should apply the matrix. i.e it is equivalent to filtering once. So, now i am thinking like i will do the filtering several times and get a lot of sub-bands i.e L0 H0 L1 H1 L2 H2 etc. i have those results in an array. (i haven't done reconstruction). Could you kindly tell me how to proceed so as to get the time and frequency detail from this. It will be a great help if i get this information. – vvv Feb 4 '16 at 5:04
• You can think of applying the matrix as half-band filtering. So If you apply the matrix to data with a sampling period $f_s$, the resulting upper band has frequency content from $\frac{f_s}{4}$ to $\frac{f_s}{2}$, and the resulting lower band has frequency content from 0 Hz to $\frac{f_s}{4}$. Discarding every second sample will reduce the sampling frequency of both of these data to $\frac{f_s}{2}$, and will alias the frequencies in the highpass band so that they are "upside down". You would continue decomposition of the lower band. So you see, DWT gives you only very low frequency resolution. – Olli Niemitalo Feb 4 '16 at 7:59
• Situation is like i am sampling a voltage signal at the rate of 500ksps. i am continuously monitoring it and if i find that a sample value has gone beyond some threshold, i will move the next (say) 1024 samples to a buffer and do a dwt on it and find the frequency and time at which those events, maybe several KHz events (if there) have occurred. Now is it possible to do this with dwt considering the frequency resolution problem you mentioned. – vvv Feb 4 '16 at 8:27