天文中存在很多干扰信号,并且很多干扰信号是远远大于背景,因为分布分散,可以考虑用中值滤波去去除这些干扰。中值滤波的原理相当于用周围的数值然后求和平均去代替每个点的数值,这样做的好处是,只要干扰比较小的情况下,可以进行有效的平滑,具体参考代码和下文的链接。
【一维】
from numpy import random
import matplotlib.pyplot as plt
from astropy.modeling import models, fitting
from astropy.stats import sigma_clip
import numpy as np
from scipy.ndimage.measurements import label
from scipy import ndimage
#自己创建一个数据
x = np.array(range(600))
y = np.random.rand(600)+10
y = y+ x*0.01
y[10:20] = y[10:20]+4
y[500:520] = y[500:520]+7 #把这两个峰当作我的干扰信号
plt.plot(x, y,label = 'before med filter')
#计算数据跟中间值的偏差,反应数据的稳定性
def getMAD(data):
med = np.ma.median(data)
mad = np.ma.median(abs(data - med))
return mad
#give data return sel_t
Nsig = 5 #设定大于5倍西伽马为干扰信号
result = ndimage.median_filter(y,size=50) #size的具体数值选取参考下文链接
diff = y - result
mad = getMAD(diff)
sel_t = abs(diff) < Nsig * mad
plt.plot(x, result, label = 'result of med filter ')
plt.plot(x[sel_t], y[sel_t], label = 'after med filter' )
plt.legend()
plt.show()
我们可以黄线是用中值滤波的结果对我们两个峰进行了有效的平滑。
【二维】
from numpy import random
import matplotlib.pyplot as plt
from astropy.modeling import models, fitting
from astropy.stats import sigma_clip
import numpy as np
from scipy.ndimage.measurements import label
from scipy import ndimage
import sys
#自己创建一个数据
array = np.random.rand(200, 300)
#自己给两个方向添加一些干扰
array[50:55, :] += 5
array[:, 70:75] += 6
#plt.imshow(array, label = 'before med filter')
#计算数据跟中间值的偏差,反应数据的稳定性
def getMAD(data):
med = np.ma.median(data)
mad = np.ma.median(abs(data - med))
return mad
#give data return sel_t
Nsig = 5 #设定大于5倍西伽马为干扰信号
result = ndimage.median_filter(array, size=100) #size的具体数值选取参考下文链接
diff = array - result
mad = getMAD(diff)
print(mad)
sel_t = abs(diff) < Nsig * mad #sel_t是中值滤波之后的数据
#plt.imshow(diff, label = 'result of med filter ')
#plt.show()
#sys.exit()
array[~sel_t] = np.nan
plt.imshow(array, label = 'after med filter' )
plt.show()
原理参考一下其他博客
https://blog.youkuaiyun.com/u013608424/article/details/72865845