当前位置:   article > 正文

图像 理想低通滤波_图像处理之滤波(下)

fft卡尔曼滤波

05eb1ea0fed6085540675fed73df1da6.png

[toc]目录

一、常规滤波

  • 低通
  • 高通
  • 带通
  • 带阻

二、非局部均值滤波

三、维纳滤波

四、卡尔曼滤波


前言

所谓滤波,其实就是从混合在一起的诸多信号中提取出所需要的信号。

信号的分类:

  • 确定型信号,可以表示为确定的时间函数,可确定其在任何时刻的量值。(具有确定的频谱);一般可通过低通、高通、带通、带阻等模拟滤波器或其他常规滤波算法实现。
  • 随机信号,不能用确定的数学关系式描述,不能预测其未来任何瞬时值,其值的变化服从统计规律。(频谱不确定,功率谱确定);根据有用信号和干扰信号的功率谱设计滤波器——维纳滤波(Wiener Filtering)或卡尔曼滤波(Kalman Filter)。

一、常规滤波

在图像处理或者计算机应用中,在正式对图像爱那个进行分析处理前一般需要一个预处理的过程。预处理是对图像作一些诸如降维、降噪的操作,主要是为后续处理提供一个体积合适、只包含所需信息的图像,通常会用到一些滤波处理手法。滤波,实际上就是信号处理,而图像本身可以看作是一个二维信号,其中像素点灰度的高低代表信号的强弱。对应的高低频的意义:

高频:图像中灰度变化强烈的点,一般是轮廓或者是噪声。

低频:图像中平坦的,灰度变化不大的点,图像中的大部各区域。

而根据图像的高频与低频的特征,可以设计相应的高通和低通滤波器,高通滤波可以检测图像中尖锐、变化明显的地方,而低通滤波可以让图像变得光滑,滤除图像中的噪声、OpenCv中提供的低通滤波有线性的均值滤波器、高斯滤波器,非线性的双边滤波器、中值滤波器;高通滤波有基于Canny,Sobel算子的各种滤波。其实很多时候低通滤波和高通滤波其实是相互矛盾的,很多时候在边缘检测前需要通过低通滤波降噪,这里就需要调节参数在保证高频的边缘不丢失的前提下尽可能多的去处图像的噪点。

这里使用频域的高通和低通滤波。

  • 低通

理想的低通滤波器的模版为:

equation?tex=H%28u%2Cv%29+%3D+%5Cbegin%7Bequation%7D+%5Cleft%5C%7B+%5Cbegin%7Barray%7D%7Blr%7D+1%2C+%5Cquad+D%28u%2C+v%29+%5Cle+D_0+%5C%5C++0%2C+%5Cquad+D%28u%2C+v%29+%5Cge+D_0+%5Cend%7Barray%7D+%5Cright.+%5Cend%7Bequation%7D

其中,

equation?tex=D_0 表示通带半径,
equation?tex=D%28u%2C+v%29 是到频谱中心的距离(欧式距离),计算公式如下:
equation?tex=D%28u%2C+v%29+%3D+%5Csqrt+%7B%28u+-+M%2F+2%29%5E2+%2B+%28v-N%2F2%29%5E2%7D ,M和N表示频谱图像的大小,
equation?tex=%28M%2F2%2C+N%2F2%29 即为频谱中心。
  1. def low_pass_filter(img, radius=100):
  2. r = radius
  3. rows, cols = img.shape
  4. center = int(rows / 2), int(cols / 2)
  5. mask = np.zeros((rows, cols, 2), np.uint8)
  6. x, y = np.ogrid[:rows, :cols]
  7. mask_area = (x - center[0]) ** 2 + (y - center[1]) ** 2 <= r * r
  8. mask[mask_area] = 1
  9. return mask

Butterworth低通滤波器为:

equation?tex=H%28u%2Cv%29+%3D+%5Cfrac%7B1%7D%7B1%2B%5B%5Cfrac%7BD%28u%2Cv%29%7D%7BD_0%7D%5D%5E%7B2n%7D%7D

fa54ce2927e67de6f636a10054a9b14d.png
  1. def Butterworth(src, d0, n, ftype):
  2. template = np.zeros(src.shape, dtype=np.float32) # 构建滤波器
  3. r, c = src.shape
  4. for i in np.arange(r):
  5. for j in np.arange(c):
  6. distance = np.sqrt((i - r/2)**2 + (j - c/2)**2)
  7. template[i, j] = 1/(1 + (distance/d0)**(2*n)) # Butterworth 滤波函数
  8. template[i, j] = np.e ** (-1 * (distance**2 / (2 * d0**2))) # Gaussian滤波函数
  9. if ftype == 'high':
  10. template = 1 - template
  11. return template

高斯低通滤波器:

equation?tex=H%28u%2Cv%29+%3D+e%5E%7B%5Cfrac%7B-D%5E2%28u%2C+v%29%7D%7B2D_0%5E2%7D%7D ,1减去低通滤波器模板即可得到高通滤波器。

e0e88eebeb638494d49c7fd975b593cf.png
  1. # 定义函数,高斯高/低通滤波模板
  2. def Gaussian(src, d0, ftype):
  3. template = np.zeros(src.shape, dtype=np.float32) # 构建滤波器
  4. r, c = src.shape
  5. for i in np.arange(r):
  6. for j in np.arange(c):
  7. distance = np.sqrt((i - r / 2) ** 2 + (j - c / 2) ** 2)
  8. template[i, j] = np.e ** (-1 * (distance ** 2 / (2 * d0 ** 2))) # Gaussian滤波函数
  9. if ftype == 'high':
  10. template = 1 - template
  11. return template
  • 带通
  • 带阻
  1. def bandreject_filters(img, r_out=300, r_in=35):
  2. rows, cols = img.shape
  3. crow, ccol = int(rows / 2), int(cols / 2)
  4. radius_out = r_out
  5. radius_in = r_in
  6. mask = np.zeros((rows, cols, 2), np.uint8)
  7. center = [crow, ccol]
  8. x, y = np.ogrid[:rows, :cols]
  9. mask_area = np.logical_and(((x - center[0]) ** 2 + (y - center[1]) ** 2 >= r_in ** 2),
  10. ((x - center[0]) ** 2 + (y - center[1]) ** 2 <= r_out ** 2))
  11. mask[mask_area] = 1
  12. mask = 1 - mask
  13. return mask

二、非局部均值滤波

图像中的像素点之间不是孤立存在的,某一点的像素与别处的像素点一定存在着某种关联,可以概括为灰度相关性和几何结构相似性。

  1. #coding:utf8
  2. import cv2
  3. import numpy as np
  4. def psnr(A, B):
  5. return 10*np.log(255*255.0/(((A.astype(np.float)-B)**2).mean()))/np.log(10)
  6. def double2uint8(I, ratio=1.0):
  7. return np.clip(np.round(I*ratio), 0, 255).astype(np.uint8)
  8. def make_kernel(f):
  9. kernel = np.zeros((2*f+1, 2*f+1))
  10. for d in range(1, f+1):
  11. kernel[f-d:f+d+1, f-d:f+d+1] += (1.0/((2*d+1)**2))
  12. return kernel/kernel.sum()
  13. def NLmeansfilter(I, h_=10, templateWindowSize=5, searchWindowSize=11):
  14. f = templateWindowSize/2
  15. t = searchWindowSize/2
  16. height, width = I.shape[:2]
  17. padLength = t+f
  18. I2 = np.pad(I, padLength, 'symmetric')
  19. kernel = make_kernel(f)
  20. h = (h_**2)
  21. I_ = I2[padLength-f:padLength+f+height, padLength-f:padLength+f+width]
  22. average = np.zeros(I.shape)
  23. sweight = np.zeros(I.shape)
  24. wmax = np.zeros(I.shape)
  25. for i in range(-t, t+1):
  26. for j in range(-t, t+1):
  27. if i==0 and j==0:
  28. continue
  29. I2_ = I2[padLength+i-f:padLength+i+f+height, padLength+j-f:padLength+j+f+width]
  30. w = np.exp(-cv2.filter2D((I2_ - I_)**2, -1, kernel)/h)[f:f+height, f:f+width]
  31. sweight += w
  32. wmax = np.maximum(wmax, w)
  33. average += (w*I2_[f:f+height, f:f+width])
  34. return (average+wmax*I)/(sweight+wmax)
  35. if __name__ == '__main__':
  36. I = cv2.imread('lena.jpg', 0)
  37. sigma = 20.0
  38. I1 = double2uint8(I + np.random.randn(*I.shape) *sigma)
  39. print u'噪声图像PSNR',psnr(I, I1)
  40. R1 = cv2.medianBlur(I1, 5)
  41. print u'中值滤波PSNR',psnr(I, R1)
  42. R2 = cv2.fastNlMeansDenoising(I1, None, sigma, 5, 11)
  43. print u'opencv的NLM算法',psnr(I, R2)
  44. R3 = double2uint8(NLmeansfilter(I1.astype(np.float), sigma, 5, 11))
  45. print u'NLM PSNR',psnr(I, R3)

尽管opencv中已经有实现,对于彩色图像,首先要先转换到CIELAB颜色空间,然后对L和AB成分分别去噪。而且据说上面的实现会比opencv自带的实现要好一些。

  1. cv2.fastNlMeansDenoising() - 使用单个灰度图像
  2. cv2.fastNlMeansDenoisingColored() - 使用彩色图像。
  3. cv2.fastNlMeansDenoisingMulti() - 用于在短时间内捕获的图像序列(灰度图像)
  4. cv2.fastNlMeansDenoisingColoredMulti() - 与上面相同,但用于彩色图像。
  5. fastNlMeansDenoisingColored( InputArray src, OutputArray dst,
  6. float h = 3, float hColor = 3,
  7. int templateWindowSize = 7, int searchWindowSize = 21)
  8. 参数:
  9. • h : 决定过滤器强度。h 值高可以很好的去除噪声,但也会把图像的细节抹去。(取 10 的效果不错)
  10. • hForColorComponents : 与 h 相同,但使用与彩色图像。(与 h 相同,10)
  11. • templateWindowSize : 奇数。(推荐值为 7)
  12. • searchWindowSize : 奇数。(推荐值为 21)

三、维纳滤波

对于运动引起的模糊,最简单的方法就是直接作逆滤波,但是逆滤波对于加性噪声特别敏感,使得恢复的图像几乎不可用。最小均方差(维纳)滤波用来去处含有噪声的模糊图像,其目标是找到未污染图像的一个估计,使得他们之间的均方误差最小,可以去除噪声,同时清晰化模糊图像。

主要参考:写的不错,有公式又证明!里面还有约束最小二乘方滤波!

https://blog.csdn.net/wsp_1138886114/article/details/95024180​blog.csdn.net
  1. import matplotlib.pyplot as plt
  2. import numpy as np
  3. from numpy import fft
  4. import math
  5. import cv2
  6. # 仿真运动模糊
  7. def motion_process(image_size, motion_angle):
  8. PSF = np.zeros(image_size)
  9. print(image_size)
  10. center_position = (image_size[0] - 1) / 2
  11. print(center_position)
  12. slope_tan = math.tan(motion_angle * math.pi / 180)
  13. slope_cot = 1 / slope_tan
  14. if slope_tan <= 1:
  15. for i in range(15):
  16. offset = round(i * slope_tan) # ((center_position-i)*slope_tan)
  17. PSF[int(center_position + offset), int(center_position - offset)] = 1
  18. return PSF / PSF.sum() # 对点扩散函数进行归一化亮度
  19. else:
  20. for i in range(15):
  21. offset = round(i * slope_cot)
  22. PSF[int(center_position - offset), int(center_position + offset)] = 1
  23. return PSF / PSF.sum()
  24. # 对图片进行运动模糊
  25. def make_blurred(input, PSF, eps):
  26. input_fft = fft.fft2(input) # 进行二维数组的傅里叶变换
  27. PSF_fft = fft.fft2(PSF) + eps
  28. blurred = fft.ifft2(input_fft * PSF_fft)
  29. blurred = np.abs(fft.fftshift(blurred))
  30. return blurred
  31. def inverse(input, PSF, eps): # 逆滤波
  32. input_fft = fft.fft2(input)
  33. PSF_fft = fft.fft2(PSF) + eps # 噪声功率,这是已知的,考虑epsilon
  34. result = fft.ifft2(input_fft / PSF_fft) # 计算F(u,v)的傅里叶反变换
  35. result = np.abs(fft.fftshift(result))
  36. return result
  37. def wiener(input, PSF, eps, K=0.01): # 维纳滤波,K=0.01
  38. input_fft = fft.fft2(input)
  39. PSF_fft = fft.fft2(PSF) + eps
  40. PSF_fft_1 = np.conj(PSF_fft) / (np.abs(PSF_fft) ** 2 + K)
  41. result = fft.ifft2(input_fft * PSF_fft_1)
  42. result = np.abs(fft.fftshift(result))
  43. return result
  44. def normal(array):
  45. array = np.where(array < 0, 0, array)
  46. array = np.where(array > 255, 255, array)
  47. array = array.astype(np.int16)
  48. return array
  49. def main(gray):
  50. channel = []
  51. img_h, img_w = gray.shape[:2]
  52. PSF = motion_process((img_h, img_w), 60) # 进行运动模糊处理
  53. blurred = np.abs(make_blurred(gray, PSF, 1e-3))
  54. result_blurred = inverse(blurred, PSF, 1e-3) # 逆滤波
  55. result_wiener = wiener(blurred, PSF, 1e-3) # 维纳滤波
  56. blurred_noisy = blurred + 0.1 * blurred.std() *
  57. np.random.standard_normal(blurred.shape) # 添加噪声,standard_normal产生随机的函数
  58. inverse_mo2no = inverse(blurred_noisy, PSF, 0.1 + 1e-3) # 对添加噪声的图像进行逆滤波
  59. wiener_mo2no = wiener(blurred_noisy, PSF, 0.1 + 1e-3) # 对添加噪声的图像进行维纳滤波
  60. channel.append((normal(blurred),normal(result_blurred),normal(result_wiener),
  61. normal(blurred_noisy),normal(inverse_mo2no),normal(wiener_mo2no)))
  62. return channel
  63. if __name__ == '__main__':
  64. image = cv2.imread('./gggg/001.png')
  65. b_gray, g_gray, r_gray = cv2.split(image.copy())
  66. Result = []
  67. for gray in [b_gray, g_gray, r_gray]:
  68. channel = main(gray)
  69. Result.append(channel)
  70. blurred = cv2.merge([Result[0][0][0], Result[1][0][0], Result[2][0][0]])
  71. result_blurred = cv2.merge([Result[0][0][1], Result[1][0][1], Result[2][0][1]])
  72. result_wiener = cv2.merge([Result[0][0][2], Result[1][0][2], Result[2][0][2]])
  73. blurred_noisy = cv2.merge([Result[0][0][3], Result[1][0][3], Result[2][0][3]])
  74. inverse_mo2no = cv2.merge([Result[0][0][4], Result[1][0][4], Result[2][0][4]])
  75. wiener_mo2no = cv2.merge([Result[0][0][5], Result[1][0][5], Result[2][0][5]])
  76. #========= 可视化 ==========
  77. plt.figure(1)
  78. plt.xlabel("Original Image")
  79. plt.imshow(np.flip(image, axis=2)) # 显示原图像
  80. plt.figure(2)
  81. plt.figure(figsize=(8, 6.5))
  82. imgNames = {"Motion blurred":blurred,
  83. "inverse deblurred":result_blurred,
  84. "wiener deblurred(k=0.01)":result_wiener,
  85. "motion & noisy blurred":blurred_noisy,
  86. "inverse_mo2no":inverse_mo2no,
  87. 'wiener_mo2no':wiener_mo2no}
  88. for i,(key,imgName) in enumerate(imgNames.items()):
  89. plt.subplot(231+i)
  90. plt.xlabel(key)
  91. plt.imshow(np.flip(imgName, axis=2))
  92. plt.show()
  93. OpenCV-Python 图像去模糊(维纳滤波,约束最小二乘方滤波)import matplotlib.pyplot as plt
  94. import numpy as np
  95. from numpy import fft
  96. import math
  97. import cv2
  98. # 仿真运动模糊
  99. def motion_process(image_size, motion_angle):
  100. PSF = np.zeros(image_size)
  101. print(image_size)
  102. center_position = (image_size[0] - 1) / 2
  103. print(center_position)
  104. slope_tan = math.tan(motion_angle * math.pi / 180)
  105. slope_cot = 1 / slope_tan
  106. if slope_tan <= 1:
  107. for i in range(15):
  108. offset = round(i * slope_tan) # ((center_position-i)*slope_tan)
  109. PSF[int(center_position + offset), int(center_position - offset)] = 1
  110. return PSF / PSF.sum() # 对点扩散函数进行归一化亮度
  111. else:
  112. for i in range(15):
  113. offset = round(i * slope_cot)
  114. PSF[int(center_position - offset), int(center_position + offset)] = 1
  115. return PSF / PSF.sum()
  116. # 对图片进行运动模糊
  117. def make_blurred(input, PSF, eps):
  118. input_fft = fft.fft2(input) # 进行二维数组的傅里叶变换
  119. PSF_fft = fft.fft2(PSF) + eps
  120. blurred = fft.ifft2(input_fft * PSF_fft)
  121. blurred = np.abs(fft.fftshift(blurred))
  122. return blurred
  123. def inverse(input, PSF, eps): # 逆滤波
  124. input_fft = fft.fft2(input)
  125. PSF_fft = fft.fft2(PSF) + eps # 噪声功率,这是已知的,考虑epsilon
  126. result = fft.ifft2(input_fft / PSF_fft) # 计算F(u,v)的傅里叶反变换
  127. result = np.abs(fft.fftshift(result))
  128. return result
  129. def wiener(input, PSF, eps, K=0.01): # 维纳滤波,K=0.01
  130. input_fft = fft.fft2(input)
  131. PSF_fft = fft.fft2(PSF) + eps
  132. PSF_fft_1 = np.conj(PSF_fft) / (np.abs(PSF_fft) ** 2 + K)
  133. result = fft.ifft2(input_fft * PSF_fft_1)
  134. result = np.abs(fft.fftshift(result))
  135. return result
  136. def normal(array):
  137. array = np.where(array < 0, 0, array)
  138. array = np.where(array > 255, 255, array)
  139. array = array.astype(np.int16)
  140. return array
  141. def main(gray):
  142. channel = []
  143. img_h, img_w = gray.shape[:2]
  144. PSF = motion_process((img_h, img_w), 60) # 进行运动模糊处理
  145. blurred = np.abs(make_blurred(gray, PSF, 1e-3))
  146. result_blurred = inverse(blurred, PSF, 1e-3) # 逆滤波
  147. result_wiener = wiener(blurred, PSF, 1e-3) # 维纳滤波
  148. blurred_noisy = blurred + 0.1 * blurred.std() *
  149. np.random.standard_normal(blurred.shape) # 添加噪声,standard_normal产生随机的函数
  150. inverse_mo2no = inverse(blurred_noisy, PSF, 0.1 + 1e-3) # 对添加噪声的图像进行逆滤波
  151. wiener_mo2no = wiener(blurred_noisy, PSF, 0.1 + 1e-3) # 对添加噪声的图像进行维纳滤波
  152. channel.append((normal(blurred),normal(result_blurred),normal(result_wiener),
  153. normal(blurred_noisy),normal(inverse_mo2no),normal(wiener_mo2no)))
  154. return channel
  155. if __name__ == '__main__':
  156. image = cv2.imread('./gggg/001.png')
  157. b_gray, g_gray, r_gray = cv2.split(image.copy())
  158. Result = []
  159. for gray in [b_gray, g_gray, r_gray]:
  160. channel = main(gray)
  161. Result.append(channel)
  162. blurred = cv2.merge([Result[0][0][0], Result[1][0][0], Result[2][0][0]])
  163. result_blurred = cv2.merge([Result[0][0][1], Result[1][0][1], Result[2][0][1]])
  164. result_wiener = cv2.merge([Result[0][0][2], Result[1][0][2], Result[2][0][2]])
  165. blurred_noisy = cv2.merge([Result[0][0][3], Result[1][0][3], Result[2][0][3]])
  166. inverse_mo2no = cv2.merge([Result[0][0][4], Result[1][0][4], Result[2][0][4]])
  167. wiener_mo2no = cv2.merge([Result[0][0][5], Result[1][0][5], Result[2][0][5]])
  168. #========= 可视化 ==========
  169. plt.figure(1)
  170. plt.xlabel("Original Image")
  171. plt.imshow(np.flip(image, axis=2)) # 显示原图像
  172. plt.figure(2)
  173. plt.figure(figsize=(8, 6.5))
  174. imgNames = {"Motion blurred":blurred,
  175. "inverse deblurred":result_blurred,
  176. "wiener deblurred(k=0.01)":result_wiener,
  177. "motion & noisy blurred":blurred_noisy,
  178. "inverse_mo2no":inverse_mo2no,
  179. 'wiener_mo2no':wiener_mo2no}
  180. for i,(key,imgName) in enumerate(imgNames.items()):
  181. plt.subplot(231+i)
  182. plt.xlabel(key)
  183. plt.imshow(np.flip(imgName, axis=2))
  184. plt.show()

四、卡尔曼滤波

Kalman filtering是一种利用线性系统状态方程,通过系统输入输出观测数据,对系统状态进行最优估计的算法。Kalman滤波在测量方差已知的情况下能够从一系列存在测量噪声的数据中,估计动态系统的状态。 A Kalman filter is an optimal estimation algorithm. (最优化自回归数据数据算法)

这都是百度百科的东西,搞得跟数学书上公式定理一样,说的都不是人话!关键是字我都认识,但我没懂什么意思。

1.1 是什么?

举个例子:在海图作业中,航海张通常以前一时刻的船位为基准,根据航向、船速和海流等一系列因素推算下一个船位,这个称之为观测船位;另外还会选择适当的方法,通过仪器得到另一个推算船位,这个称之为推算船位;观测船位和推算船位一般都不重合,航海长需要通过分析和判断选择一个可靠的船位,作为船舰当前的位置。

由此得出卡尔曼滤波思想:以

equation?tex=K-1 时刻的最有估计
equation?tex=x_%7Bk-1%7D 为准,预测
equation?tex=K 时刻的状态变量
equation?tex=%5Chat+x_%7Bk%2Fk-1%7D ,同时又对该状态进行观测,得到观测变量
equation?tex=z_k ,再在预测与预测之间进行分析,或者说是以观测量对预测量进行修正,从而得到
equation?tex=K 时刻的最有状态估计
equation?tex=x_k

卡尔曼滤波算法在控制领域有着极广泛的应用,比如自动驾驶汽车,想象一下,一个雷达传感器告诉你另一辆车距离15米,一个激光传感器说车辆距离20米 ,那么如何协调这些传感器测量。再比如,在发动机燃油喷射控制中,可以应用扩展的卡尔曼滤波理论研究瞬态工况下发动机循环进气量的最有估计算法;在雷达中,人们最感兴趣的是跟踪目标,但目标的位置、速度、加速度的测量往往在任何时候都有噪声,而卡尔曼滤波则是利用目标的动态信息,设法去除噪声,得到一个关于目标位置的最好估计。

1.2 怎么回事?

再举个例子:假设要研究一个房间的温度,以一分钟位时间单位,根据经验判断,这个房间的温度是恒定的,但是对于我们的经验不是完全相信,可能存在上下几度的偏差,我们把该偏差看作是高斯白噪声。另外,在房间里放一个温度计,而温度计也不准确,测量值会与实际值存在偏差,我们也把这偏差看作是高斯白噪声。那么现在,我们要根据经验温度和温度计的测量值以及他们各自的噪声来估算出放房间的实际温度。

那么接下来该如何解决呢? 假如我们要估算

equation?tex=k 时刻的实际温度,首先要根据
equation?tex=k-1 时刻的温度值,来预测
equation?tex=k 时刻的温度(
equation?tex=K 时刻的经验温度)。因为我们的经验认为温度是恒定的,所以会得到
equation?tex=k 时刻的温度和
equation?tex=k-1 时刻是一样的,假设是23度,同时该值(预测值)的高斯白噪声为5度,(5是这样得到的。如果k-1时刻估算出的最优温度值的偏差是3,对自己预测的不确定度是4度,他们平方相加再开方,就是5),而温度计得到的温度值为25度,同时该值的偏差为4,此时,对于K时刻房间的温度值有两个:估计值23度和测量值25度,那么究竟相信谁?用均方根误差判断,
equation?tex=H%5E2+%3D+%5Cfrac%7B5%5E2%7D%7B5%5E2+%2B+4%5E2%7D+%5CRightarrow+H+%3D+0.78

所以估计K时刻最优温度值为

equation?tex=23%2B0.78%2A%2825-23%29%3D24.56 度,得到了K时刻的最优温度,下一步就是对K+1时刻的温度值进行最优估算,需要K时刻的最优温度的偏差,
equation?tex=%5Csqrt+%7B%281-H%29+%2A5%5E2%7D+%3D+2.35

fd8a67d6e9b2e5b6eb4c68c0b54fbd9c.png

e7cd62d2b0759bf24353fc2e85bd68be.png
无控制离散型卡尔曼滤波基本方程

ccf57598e97f0ea3494a680b861b0811.png

112f38093ea4eb43df1adc057297a630.png
带有控制的离散型卡尔曼滤波基本方程

详情参照:[卡尔曼滤波器分类及其基本公式](卡尔曼滤波器分类及基本公式 - 百度文库)

就这样,卡尔曼滤波就能不断吧均方误差递归,从而估算出最优的温度值,运行速度快,且只保留上一时刻的协方差。

总而言之,Kalman滤波用在当测量值与模型预测值均不准确的情况下,用来计算预测真值的一种滤波算法,在目标识别与追踪任务中经常用到。

1.3 python代码实现

  1. def KalmanFilter(z, n_iter=20):
  2. # 这里是假设A=1,H=1的情况
  3. # intial parameters
  4. sz = (n_iter,) # size of array
  5. # Q = 1e-5 # process variance
  6. Q = 1e-6 # process variance
  7. # allocate space for arrays
  8. xhat = numpy.zeros(sz) # a posteri estimate of x
  9. P = numpy.zeros(sz) # a posteri error estimate
  10. xhatminus = numpy.zeros(sz) # a priori estimate of x
  11. Pminus = numpy.zeros(sz) # a priori error estimate
  12. K = numpy.zeros(sz) # gain or blending factor
  13. R = 0.1 ** 2 # estimate of measurement variance, change to see effect
  14. # intial guesses
  15. xhat[0] = 0.0
  16. P[0] = 1.0
  17. A = 1
  18. H = 1
  19. for k in range(1, n_iter):
  20. # time update
  21. xhatminus[k] = A * xhat[k - 1] # X(k|k-1) = AX(k-1|k-1) + BU(k) + W(k),A=1,BU(k) = 0
  22. Pminus[k] = A * P[k - 1] + Q # P(k|k-1) = AP(k-1|k-1)A' + Q(k) ,A=1
  23. # measurement update
  24. K[k] = Pminus[k] / (Pminus[k] + R) # Kg(k)=P(k|k-1)H'/[HP(k|k-1)H' + R],H=1
  25. xhat[k] = xhatminus[k] + K[k] * (z[k] - H * xhatminus[k]) # X(k|k) = X(k|k-1) + Kg(k)[Z(k) - HX(k|k-1)], H=1
  26. P[k] = (1 - K[k] * H) * Pminus[k] # P(k|k) = (1 - Kg(k)H)P(k|k-1), H=1
  27. return xhat

每日一句毒鸡汤:

牛逼的算法往往都是来源一个很简单的思想所演化出来的!

继续坚持!加油~

声明:本文内容由网友自发贡献,不代表【wpsshop博客】立场,版权归原作者所有,本站不承担相应法律责任。如您发现有侵权的内容,请联系我们。转载请注明出处:https://www.wpsshop.cn/w/很楠不爱3/article/detail/117282
推荐阅读
  

闽ICP备14008679号