Previous Up Next

19.4.5  Discrete wavelet transform

The dwt and idwt commands are used for computing the discrete wavelet transform (DWT) of a signal and for computing the inverse transformation, i.e. reconstructing the original signal from the transform. These commands require the GSL library.

Examples

In the first example (adapted from the GSL documentation, see https://www.gnu.org/software/gsl/doc/html/dwt.html), we load a sample of size 256 from the MIT-BIH Arrhythmia Database (see https://physionet.org/content/mitdb/1.0.0/) and attempt to remove the noisy information.

data:=mid(col(csv2gen("/home/luka/Documents/MIT-BIH.csv","\t"),2),300,256):; listplot(data)

Now we transform the data using dwt and select the 20 largest components. We set other elements to zero and transform the result back using idwt.

tdata:=dwt(data):; p:=reverse(sortperm(abs(tdata))):; for k from 20 to 255 do tdata[p[k]]:=0; od:; res:=idwt(tdata):; listplot(res)

In the second example, we attempt to reduce the noise in a noisy image.

img:=image("/home/luka/Pictures/noisy.jpg")(grey)
     
an image of size 300×300 (grayscale)           

We transform the image matrix and replace each component smaller than 55 (by absolute value) in the resulting matrix to zero, before transforming it back to obtain a noise-reduced image. (The threshold value was obtained by trial and error.) Note that dwt enlarges the input matrix to the size 512× 512, so we crop the output of idwt back to 300× 300. We also make sure that no component gets below zero or above 255.

tdata:=dwt(img[0],image):; sdata:=threshold(w,55.0=0,'<=',abs=true):; nr:=image(1,threshold(round(subMat(idwt(sdata,image),0,0,299,299)),[0,255])):;

Finally, we display the original image and the processed image next to each other for comparison.

display(img,0); display(nr,320); legend(-20i,"original"); legend(320-20i,"denoised")

Previous Up Next