Estimating background with SWT

[Lech Piotrowski]

Hello,

I am looking for significant signal peaks on top of a slowly variable background. An example of the data that I am analysing looks like that:

In principle to get the best SNR I would:

1. Subtract the background

2. Remove the point-to-point variability by thresholding high levels of the wavelet decomposition.

While (2) seems to work reasonably, I have troubles estimating the background. I thought that in principle when I do:

s_swt = pywt.swt(s, "sym3", level=6)

and then reconstruct the signal zeroing the first approximates s_swt[0][0] I should get the signal on a "flat" background. However, it

The background seems to be more "bent" towards the ends, not to mention a complete change of scale.

I guess this approach is somehow completely wrong. I would be grateful if someone could tell me why and what would be the proper approach.

I also tried to use non-reconstructed first approximates as background. The shape more or less fits, but not the scale, even though here I've used "norm=True" for the decomposition. Here I plot it as a line on top of the data:

On the other hand, signal reconstructed just from s_swt[0][0] (all other details/approximates zeroed) do not match the background at all, which is most puzzling for me:

I would appreciate any help.

[Gregory Lee]

Hi Lech,

I think you are generally on the right track. At least for the bundled ecg data in pywt.data.ecg(), it does look like it helps with background removal. Try for example:

```Python

import numpy as np

import pywt

x = pywt.data.ecg()

c = pywt.swt(x, 'sym3', level=6)

r = pywt.iswt(c, 'sym3')

```

If you plot r vs. x in the above, it does look like the background has become flatter since low frequency components were removed.

If you are doing a round trip transform like above, it should not matter which normalization option you choose as long as you use the same one for BOTH the `swt` and `iswt`, the overall scaling should be consistent. It might be that for the relatively short data you used, boundary effects are coming into play or you may need to adjust the number of wavelet levels in the decomposition/reconstruction.

A second approach you may want to try is using skimage.restoration.rolling_ball which is in scikit-image >= 0.18 (see: https://scikit-image.org/docs/stable/auto_examples/segmentation/plot_rolling_ball.html). Despite being designed for images, the implementation is n-dimensional and works for 1D data. Here is an example for the same data above:

```Python

from skimage import restoration

from matplotlib import pyplot as plt

import pywt

x = pywt.data.ecg()

background = restoration.rolling_ball(x, radius=80)

plt.figure()

plt.plot(x)

plt.plot(x - background)

# a smaller radius would remove some of the smoother peaks in the ECG as well

background2 = restoration.rolling_ball(x, radius=10)

plt.plot(x - background2)

plt.legend(('original', 'radius 80', 'radius 10'))

```

Here, the radius can be used to control the size of features that get considered part of the background. You can also try the peak_local_max function in scikit-image to detect peaks (https://scikit-image.org/docs/dev/api/skimage.feature.html#skimage.feature.peak_local_max).

Cheers,

Greg

--
You received this message because you are subscribed to the Google Groups "PyWavelets" group.
To unsubscribe from this group and stop receiving emails from it, send an email to pywavelets+...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/pywavelets/80238a05-ef25-41d3-82fa-78ced91f69d5n%40googlegroups.com.

[Lech Piotrowski]

Thank you very much for your answer.

However, I see some similar effects on the ecg that I could see on my pictures:

1. The overall trend at the beginning (up to ~200) is removed, but some trends are added later. The spectrum around 400 and 700 develops a bump after filtering and was much flatter before filtering.

2. The whole scale is shifted. That may actually be the expected behaviour.

So, point 1. is my main worry. Perhaps it is as it should be with wavelets, but honestly I am not sure I understand why those bumps develop. Perhaps it is an effect of convolution of a very wide wavelet with those high peaks...

And thanks for the rolling ball. I'll definitely take a look.