0
$\begingroup$

Main Problem: How can you inverse Wavelet Transforms

(Using the data given by signal.scipy.cwt)

I was wondering if anyone understands the scipy.signal.cwt() function well enough to use it to remove backgrounds and noise from data. I have seen where Matlab has an inverse continuous wavelet transform function which will return the original form of the data by inputting the wavelet transform, although you can filter out the slices you don't want.

MATALAB inverse cwt funciton

Since scipy doesn't appear to have the same function, I have been trying to figure out how to get the data back in the same form, while removing the noise and background. How do I do this? I noticed that the wavelet has negative values, when the original data is made of Gaussian-like feautures, so I tried squaring it, but this gives me values way to large and not quite right. I tried reading about the transforms but everything is very complicated which was difficult to parse in order to figure out what I needed to do.

Here is what I have been trying:

# Compute the wavelet transform
# I can't figure out what the width is or does?
widths = range(1,11) 
# Ricker is 2nd derivative of Gaussian
# (*close* to what *most* of the features are in my data)
# (They're actually Lorentzians and Breit-Wigner-Fano lines)
cwtmatr = signal.cwt(xy['y'], signal.ricker, widths)

# Maybe we multiple by the original data? and square?
WT_to_original_data = (xy['y'] * cwtmatr)**2

And here is a fully compilable short script to show you the type of data I am trying to get and what I have etc.:

import numpy as np
from scipy import signal
import matplotlib.pyplot as plt

# Make some random data with peaks and noise
def make_peaks(x):
    bkg_peaks = np.array(np.zeros(len(x)))
    desired_peaks = np.array(np.zeros(len(x)))
    # Make peaks which contain the data desired
    # (Mid range/frequency peaks)
    for i in range(0,10):
        center = x[-1] * np.random.random() - x[0]
        amp = 60 * np.random.random() + 10
        width = 10 * np.random.random() + 5
        desired_peaks += amp * np.e**(-(x-center)**2/(2*width**2))
    # Also make background peaks (not desired)
    for i in range(0,3):
        center = x[-1] * np.random.random() - x[0]
        amp = 40 * np.random.random() + 10
        width = 100 * np.random.random() + 100
        bkg_peaks += amp * np.e**(-(x-center)**2/(2*width**2))
    return bkg_peaks, desired_peaks

x = np.array(range(0, 1000))
bkg_peaks, desired_peaks = make_peaks(x)
y_noise = np.random.normal(loc=30, scale=10, size=len(x))
y = bkg_peaks + desired_peaks + y_noise
xy = np.array( zip(x,y), dtype=[('x',float), ('y',float)])

# Compute the wavelet transform
# I can't figure out what the width is or does?
widths = range(1,11) 
# Ricker is 2nd derivative of Gaussian
# (*close* to what *most* of the features are in my data)
# (They're actually Lorentzians and Breit-Wigner-Fano lines)
cwtmatr = signal.cwt(xy['y'], signal.ricker, widths)

# Maybe we multiple by the original data? and square?
WT = (xy['y'] * cwtmatr)**2


# plot the data and results
fig = plt.figure()
ax_raw_data = fig.add_subplot(4,3,1)
ax = {}
for i in range(0, 11):
    ax[i] = fig.add_subplot(4,3, i+2)

ax_desired_transformed_data = fig.add_subplot(4,3,12)

ax_raw_data.plot(xy['x'], xy['y'], 'g-')  
for i in range(0,10):
    ax[i].plot(xy['x'], WT[i])

ax_desired_transformed_data.plot(xy['x'], desired_peaks, 'k-')

fig.tight_layout()
plt.show()

This script will output this image:

output

Where the first plot is the raw data, the middle plots are the wavelet transforms and the last plot is what I want to get out as the processed (background and noise removed) data.

Does anyone have any suggestions? Thank you so much for the help.

$\endgroup$
  • $\begingroup$ The question is interesting, but you should really build a solid background before attempting to implement such an algorithm. If you want to kick-start the things - read about wavelet thresholding. $\endgroup$ – Sektor Jan 28 '14 at 17:53
0
$\begingroup$

I ended up finding a package which provides an inverse wavelet transform function called mlpy. The function is mlpy.wavelet.uwt. This is the snippet which may interest people if they are trying to do noise or background removal:

# Make 2**n amount of data
new_y, bool_y = wave.pad(xy['y'])
orig_mask = np.where(bool_y==True)

# wavelet transform parameters
levels = 8
wf = 'h'
k = 2

# Remove Noise first
# Wave transform
wt = wave.uwt(new_y, wf, k, levels)
# Matrix of the difference between each wavelet level and the original data
diff_array = np.array([(wave.iuwt(wt[i:i+1], wf, k)-new_y) for i in range(len(wt))])
# Index of the level which is most similar to original data (to obtain smoothed data)
indx = np.argmin(np.sum(diff_array**2, axis=1))
# Use the wavelet levels around this region
noise_wt = wt[indx:indx+1]
# smoothed data in 2^n length
new_y = wave.iuwt(noise_wt, wf, k)

# Background Removal
error = 10000
errdiff = 100
i = -1
iter_y_dict = {0:np.copy(new_y)}
bkg_approx_dict = {0:np.array([])}
while abs(errdiff)>=1*10**-24:
    i += 1
    # Wave transform
    wt = wave.uwt(iter_y_dict[i], wf, k, levels)

    # Assume last slice is lowest frequency (background approximation)
    bkg_wt = wt[-3:-1]
    bkg_approx_dict[i] = wave.iuwt(bkg_wt, wf, k)

    # Get the error
    errdiff = error - sum(iter_y_dict[i] - bkg_approx_dict[i])**2
    error = sum(iter_y_dict[i] - bkg_approx_dict[i])**2

    # Make every peak higher than bkg_wt
    diff = (new_y - bkg_approx_dict[i])
    peak_idxs_to_remove = np.where(diff>0.)[0]
    iter_y_dict[i+1] = np.copy(new_y)
    iter_y_dict[i+1][peak_idxs_to_remove] = np.copy(bkg_approx_dict[i])[peak_idxs_to_remove]

# new data without noise and background
new_y = new_y[orig_mask]
bkg_approx = bkg_approx_dict[len(bkg_approx_dict.keys())-1][orig_mask]
new_data = diff[orig_mask] 

Here a complete compilable script which can produce output:

import numpy as np
from scipy import signal
import matplotlib.pyplot as plt
import mlpy.wavelet as wave

# Make some random data with peaks and noise
############################################################
def gen_data():
    def make_peaks(x):
        bkg_peaks = np.array(np.zeros(len(x)))
        desired_peaks = np.array(np.zeros(len(x)))
        # Make peaks which contain the data desired
        # (Mid range/frequency peaks)
        for i in range(0,10):
            center = x[-1] * np.random.random() - x[0]
            amp = 100 * np.random.random() + 10
            width = 10 * np.random.random() + 5
            desired_peaks += amp * np.e**(-(x-center)**2/(2*width**2))
        # Also make background peaks (not desired)
        for i in range(0,3):
            center = x[-1] * np.random.random() - x[0]
            amp = 80 * np.random.random() + 10
            width = 100 * np.random.random() + 100
            bkg_peaks += amp * np.e**(-(x-center)**2/(2*width**2))
        return bkg_peaks, desired_peaks

    # make x axis
    x = np.array(range(0, 1000))
    bkg_peaks, desired_peaks = make_peaks(x)
    avg_noise_level = 30
    std_dev_noise = 10
    size = len(x)
    scattering_noise_amp = 100
    scat_center = 100
    scat_width = 15
    scat_std_dev_noise = 100
    y_scattering_noise = np.random.normal(scattering_noise_amp, scat_std_dev_noise, size) * np.e**(-(x-scat_center)**2/(2*scat_width**2))
    y_noise = np.random.normal(avg_noise_level, std_dev_noise, size) + y_scattering_noise
    y = bkg_peaks + desired_peaks + y_noise
    xy = np.array( zip(x,y), dtype=[('x',float), ('y',float)])
    return xy
# Random data Generated
#############################################################

xy = gen_data()

# Make 2**n amount of data
new_y, bool_y = wave.pad(xy['y'])
orig_mask = np.where(bool_y==True)

# wavelet transform parameters
levels = 8
wf = 'h'
k = 2

# Remove Noise first
# Wave transform
wt = wave.uwt(new_y, wf, k, levels)
# Matrix of the difference between each wavelet level and the original data
diff_array = np.array([(wave.iuwt(wt[i:i+1], wf, k)-new_y) for i in range(len(wt))])
# Index of the level which is most similar to original data (to obtain smoothed data)
indx = np.argmin(np.sum(diff_array**2, axis=1))
# Use the wavelet levels around this region
noise_wt = wt[indx:indx+1]
# smoothed data in 2^n length
new_y = wave.iuwt(noise_wt, wf, k)

# Background Removal
error = 10000
errdiff = 100
i = -1
iter_y_dict = {0:np.copy(new_y)}
bkg_approx_dict = {0:np.array([])}
while abs(errdiff)>=1*10**-24:
    i += 1
    # Wave transform
    wt = wave.uwt(iter_y_dict[i], wf, k, levels)

    # Assume last slice is lowest frequency (background approximation)
    bkg_wt = wt[-3:-1]
    bkg_approx_dict[i] = wave.iuwt(bkg_wt, wf, k)

    # Get the error
    errdiff = error - sum(iter_y_dict[i] - bkg_approx_dict[i])**2
    error = sum(iter_y_dict[i] - bkg_approx_dict[i])**2

    # Make every peak higher than bkg_wt
    diff = (new_y - bkg_approx_dict[i])
    peak_idxs_to_remove = np.where(diff>0.)[0]
    iter_y_dict[i+1] = np.copy(new_y)
    iter_y_dict[i+1][peak_idxs_to_remove] = np.copy(bkg_approx_dict[i])[peak_idxs_to_remove]

# new data without noise and background
new_y = new_y[orig_mask]
bkg_approx = bkg_approx_dict[len(bkg_approx_dict.keys())-1][orig_mask]
new_data = diff[orig_mask] 

##############################################################
# plot the data and results
fig = plt.figure()

ax_raw_data = fig.add_subplot(121)
ax_WT = fig.add_subplot(122)

ax_raw_data.plot(xy['x'], xy['y'], 'g')
for bkg in bkg_approx_dict.values():
    ax_raw_data.plot(xy['x'], bkg[orig_mask], 'k')

ax_WT.plot(xy['x'], new_data, 'y')


fig.tight_layout()
plt.show()

And here is the output I am getting now: Shifting Background Removal As you can see, there is still a problem with the background removal (it shifts to the right after each iteration), but it is a different question which I will address here.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.