Signal Denoising with Pywavelets

[DavidG]

Can any one provide an example of signal denoiseing using the
functions from PyWavelets?
I'm not sure how to proceed. Here's my attempt:

# ******************Code starts here*************************
import numpy as np
import pywt
import random
import pylab as plt


def main():

(cA, cD) = pywt.dwt([1,2,3,4,5,6,7,8,9,10], 'db2', 'sp1')
for i,coeff in enumerate(cA):
print i,coeff,cD[i]
print pywt.idwt(cA, cD, 'db2', 'sp1')
#[ 1. 2. 3. 4. 5. 6.]

# construct a sin wave and add noise
puresig=np.sin(np.arange(0,60,.1))
noisesig=np.random.rand(len(puresig))-.5 # subtract .5 offsets the
noise values to the interval -.5 to .5
testsig=puresig+.2*noisesig # add the noise to the signal
print 'length of signal vectors
=',len(puresig),len(noisesig),len(testsig)

# If I understood correctly this should return only the first 3
coefficients of the wavelet transform for the test signal
coeffs = pywt.wavedec(testsig, 'db2', level=2)
# the following should return the reconstructed "denoised" signal
recon2= pywt.waverec(coeffs, 'db2')
print "db2 ",len(coeffs)
print

# when I plot the data its still very noisy is this correct????
plt.plot(puresig,'k')
plt.plot(testsig,'b')

plt.plot(recon2,'g')

plt.show()

if __name__ == '__main__':
main()


# ******************Code ends here*************************

Any help or advice would be appreciated.

[Michael Marino]

Hi David-

Check out the attached code, it should give you an idea of how to do this. The iswt function is the inverse swt function which is (I don't think yet?) in pyWavelets but probably necessary for what you want to do. The dwt will introduce Gibbs artifacts into your data, so it's preferential to use a time-invariant wavelet transformation. (Donoho has published several articles on this which could be of some interest.) The actual thresholding used is up to you and should be dependent upon your particular application. I would encourage you to look through the many research articles available which go through this to avoid applying the following code as a black box.

Begin code ********************

def iswt(coefficients, wavelet):
"""
M. G. Marino to complement pyWavelets' swt.
Input parameters:

coefficients
approx and detail coefficients, arranged in level value
exactly as output from swt:
e.g. [(cA1, cD1), (cA2, cD2), ..., (cAn, cDn)]

wavelet
Either the name of a wavelet or a Wavelet object

"""
output = coefficients[0][0].copy() # Avoid modification of input data

#num_levels, equivalent to the decomposition level, n
num_levels = len(coefficients)
for j in range(num_levels,0,-1):
step_size = int(math.pow(2, j-1))
last_index = step_size
_, cD = coefficients[num_levels - j]
for first in range(last_index): # 0 to last_index - 1

# Getting the indices that we will transform
indices = numpy.arange(first, len(cD), step_size)

# select the even indices
even_indices = indices[0::2]
# select the odd indices
odd_indices = indices[1::2]

# perform the inverse dwt on the selected indices,
# making sure to use periodic boundary conditions
x1 = pywt.idwt(output[even_indices], cD[even_indices], wavelet, 'per')
x2 = pywt.idwt(output[odd_indices], cD[odd_indices], wavelet, 'per')

# perform a circular shift right
x2 = numpy.roll(x2, 1)

# average and insert into the correct indices
output[indices] = (x1 + x2)/2.

return output

def apply_threshold(output, scaler = 1., input=None):

"""
output
approx and detail coefficients, arranged in level value
exactly as output from swt:
e.g. [(cA1, cD1), (cA2, cD2), ..., (cAn, cDn)]
scaler
float to allow runtime tuning of thresholding
input
vector with length len(output). If not None, these values are used for thresholding
if None, then the vector applies a calculation to estimate the proper thresholding
given this waveform.
"""

for j in range(len(output)):
cA, cD = output[j]
if input is None:
dev = numpy.median(numpy.abs(cD - numpy.median(cD)))/0.6745
thresh = math.sqrt(2*math.log(len(cD)))*dev*scaler
else: thresh = scaler*input[j]
cD = pywt.thresholding.hard(cD, thresh)
output[j] = (cA, cD)

def measure_threshold(output, scaler = 1.):
"""
output
approx and detail coefficients, arranged in level value
exactly as output from swt:
e.g. [(cA1, cD1), (cA2, cD2), ..., (cAn, cDn)]
scaler
float to allow runtime tuning of thresholding

returns vector of length len(output) with treshold values

"""
measure = []
for j in range(len(output)):
cA, cD = output[j]
dev = numpy.median(numpy.abs(cD - numpy.median(cD)))/0.6745
thresh = math.sqrt(2*math.log(len(cD)))*dev*scaler
measure.append(thresh)
return measure

End code ********************

Proper use should look like:

output = pywt.swt(vec, wl_trans, level=level)
apply_threshold(output, 0.8, threshold_list)
final = iswt(output, wl_trans)

where threshold_list might be the return of measure_threshold, and wl_trans is either a string or a Wavelet.

In general, I would suggest estimating the thresholds using some large ensemble of noise pulses.

Cheers,
Mike

> --
> You received this message because you are subscribed to the Google Groups "PyWavelets" group.
> To post to this group, send email to pywav...@googlegroups.com.
> To unsubscribe from this group, send email to pywavelets+...@googlegroups.com.
> For more options, visit this group at http://groups.google.com/group/pywavelets?hl=en.
>

[Rick Muller]

Apologies both newbie questions and for replying to a 6 month-old thread, but I'm trying to do something very similar to the OPs intent, and thought maybe it was more useful to add to the thread.

I'm working with a group of experimental physicists, trying to help pull a signal out of the noise in an electronic device test.

What's the difference between using SWT and wavedec here? My own first guess at how to denoise the signal was to perform a multilevel decomposition, and then use some thresholding function at each level to remove the noise. However, looking at the wikipedia page for the SWT, it seems like there is something multilevel going on there as well. Can someone help explain the basics to me?