移位不变小波变换的移位

信息处理 小波 Python 转换 阈值
2022-02-08 01:07:09

主要问题:为什么迭代小波/逆小波变换会导致未抽取(移位不变)小波变换沿 x 轴移位?

我正在尝试使用类似于我在一篇文章中找到的这种方法的迭代小波变换方法从信号中删除背景:

文章描述

但是,我从下面描述的 python 程序中收到了这个输出:

换档

我不明白为什么逆小波将每次迭代都向右移动。什么可能导致这种情况?

这是我用来生成此输出的脚本:

import numpy as np
import matplotlib.pyplot as plt
import mlpy.wavelet as wave

# This fucntion should be fine 
# 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
#############################################################


#############################################################
# Wavelet Transformations
#############################################################

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] 


#############################################################
# This part should be fine
# 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()
0个回答
没有发现任何回复~