-
Notifications
You must be signed in to change notification settings - Fork 0
/
db_esti_bg.py
91 lines (74 loc) · 2.55 KB
/
db_esti_bg.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
# -*- coding: utf-8 -*-
"""
Created on Mon Sep 9 10:44:03 2019
@author: xiaozhang
"""
## estimation_background from falcon2d estimation_background.m
from pywt import wavedec,waverec
import numpy as np
import matplotlib.pyplot as plt
t = np.linspace(-100,100,50)
bg = 5
y = 50*np.exp(-t**2/(2*2^2)) + bg
plt.plot(t,y)
iter = 5
x_filt = y
th = 1
for i in range(iter):
#wavelet transform
coeffs = wavedec(x_filt,'db3',level=3)
cA = coeffs[0]
cD = []
# 令高频信息为0,也就是令cDn全部设置为0
for j in range(1,len(coeffs)):
a = np.squeeze( np.array([np.zeros([np.size(coeffs[j]),1 ])]) )
cD.append(a)
cD.insert(0,cA)
# inverse wavelet transform by only using low-freq components
x_new = waverec(cD, 'db3')
if th > 0:
# cut off values over current estimated background level.
eps = np.sqrt(abs(x_filt))/2
ind = y > (x_new + eps)
x_filt[ind] = x_new[ind]+eps[ind]
# re-estimate background
coeffs = wavedec(x_filt,'db3',level=3)
cA = coeffs[0]
cD = []
for j in range(1,len(coeffs)): #排除cA5,cDn全部设置为0
a = np.squeeze( np.array([np.zeros([np.size(coeffs[j]),1 ])]) )
cD.append(a)
cD.insert(0,cA)
x_new = waverec(cD, 'db3')
def est_bg(y, iter):
x_filt = y
th = 1
for i in range(iter):
#wavelet transform
coeffs = wavedec(x_filt,'db6',level=5)
cA = coeffs[0]
cD = []
# 令高频信息为0,也就是令cDn全部设置为0
for j in range(1,len(coeffs)):
a = np.squeeze( np.array([np.zeros([np.size(coeffs[j]),1 ])]) )
cD.append(a)
cD.insert(0,cA)
# inverse wavelet transform by only using low-freq components
x_new = waverec(cD, 'db6')
if th>0:
#cut off values over current estimated background level.
eps = np.sqrt(abs(x_filt))/2
ind = y > (x_new + eps)
x_filt[ind] = x_new[ind]+eps[ind]
# re-estimate background
coeffs = wavedec(x_filt,'db6',level=5)
cA = coeffs[0]
cD = []
# 令高频信息为0,也就是令cDn全部设置为0
for j in range(1,len(coeffs)):
a = np.squeeze( np.array([np.zeros([np.size(coeffs[j]),1 ])]) )
cD.append(a)
cD.insert(0,cA)
x_new = waverec(cD, 'db6')
return x_new
#est_bg(y,5)