ndimage.median_filter中值滤去除干扰信号

天文中存在很多干扰信号,并且很多干扰信号是远远大于背景,因为分布分散,可以考虑用中值滤波去去除这些干扰。中值滤波的原理相当于用周围的数值然后求和平均去代替每个点的数值,这样做的好处是,只要干扰比较小的情况下,可以进行有效的平滑,具体参考代码和下文的链接。

【一维】

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

评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值