信息处理 - 估计被高频噪声破坏的低频信号 - 吾爱随笔录

估计被高频噪声破坏的低频信号

信息处理离散信号噪音估计

2022-01-27 02:24:05

这是对如何为发表研究工作选择期刊？问题。假设我们对被高频噪声破坏的低频信号进行离散采样。我们知道信号的频率没有高于奈奎斯特速率的 20%。我们还知道所有噪声都在奈奎斯特速率的 50% 或以上，并且低于奈奎斯特速率。我将使用以下作为具体示例进行实验。

\begin{aligned} measurement (t) = signal (t) + noise (t) \\ signal (t) = 4.4 + \sin (0.06 π t - 0.1) + \sin (0.14 π t + 0.07) + \sin (0.2 π t + 0.4) \\ noise (t) = \sin (0.06 + 0.5 π t) + \sin (0.1 - 0.68 π t) - \sin (0.04 - 0.74 π t) + \sin (0.03 - 0.93 π t) \end{aligned}

$\begin{align}&\textrm{measurement}(t) = \textrm{signal}(t) + \textrm{noise}(t)\\ &\textrm{signal}(t) = 4.4+\sin(0.06\pi t-0.1) + \sin(0.14\pi t+0.07) + \sin(0.2\pi t+0.4)\\ &\textrm{noise}(t) = \sin(0.06+0.5\pi t) + \sin(0.1-0.68\pi t) - \sin(0.04-0.74\pi t) + \sin(0.03-0.93\pi t) \end{align}$

)的整数值进行采样。下面是测量样本和连续信号的图。 $t$

使用当前和过去的测量 )的整数值处估计我能够使用下面的卡尔曼滤波器对其进行过滤。 $\textrm{signal}(t)$ $t$

信号和卡尔曼滤波器输出的图在下一个图中。最后的图显示了卡尔曼滤波器输出中的误差，其最大值为 1.865，RMS 误差为 0.7895。我们可以制作一个比上面的卡尔曼滤波器做得更好的卡尔曼滤波器吗？什么已知的过滤方法会给出最好的估计，估计会有多好？

25 年前我了解了离散低通滤波器，我一直认为上述问题是它们的常见应用。当时不知道这个问题叫估计，也不知道卡尔曼滤波器。然后我最近拿了这个当然，佐治亚理工学院提供的另一门课程详细介绍了估计算法。然而，佐治亚理工学院的课程没有提到使用离散低通滤波器进行估计。我查了几本关于估计的书，没有发现任何关于离散低通滤波器的信息。我和一位从事数据融合研究 20 年的计算机科学家一起工作，他对 z 变换、传递函数等一无所知。为什么估计界忽略了使用离散低通滤波器进行估计？据我所知，这是估计上述信号的最佳方法。

1个回答

要回答你的最后一个问题：

为什么估计界忽略了使用离散低通滤波器进行估计？据我所知，这是估计上述信号的最佳方法。

那是因为你给他们提供了错误的信号模型。

TL;DR：为工作使用正确的工具！

血腥细节

任何时候你从卡尔曼滤波器开始，你都假设信号（measurement在上面的方程中）是这样生成的：其中和和是零均值，高斯分布，独立具有协方差的随机变量至少对于你引用的方程。KF 的更一般应用是可能的，但通常不会这样做。

x [k + 1] = A x [k] + u [k] measurement [k] = H x [k] + v [k]

$x[k+1] = \mathbf{A} x[k] + u[k]\\ \mbox{measurement}[k] = \mathbf{H} x[k] + v[k]$

A = [\begin{array}{cc} 1 & 1 \\ 0 & 1 \end{array}] H = [\begin{array}{cc} 1 & 0 \end{array}]

$\mathbf{A} = \left [ \begin{array}{cc} 1 & 1 \\ 0 & 1 \end{array} \right]\\ \mathbf{H} = \left [ \begin{array}{cc} 1 & 0 \end{array} \right]\\$

u

$u$

v

$v$

Q = [\begin{array}{cc} 1.6 & 2.4 \\ 2.4 & 4.8 \end{array}] \times 10^{- 9} R = 2.35 \times 10^{- 6}

$\mathbf{Q} = \left [ \begin{array}{cc} 1.6 & 2.4 \\ 2.4 & 4.8 \end{array} \right] \times 10^{-9}\\ \mathbf{R} = 2.35 \times 10^{-6}\\$

所以，可以做的一件事来改进你的估计是确保和更接近现实。从外观上看，您的过程噪声（ ] ）具有大致相同的方差。让他们更亲近。 $\mathbf{Q}$ $\mathbf{R}$ $u[k]$ $v[k]$

更重要的是，您的信号与您的噪声相关，因为您有谐波噪声，而不是随机噪声。更糟糕的是：您有同步谐波噪声（意味着噪声与您的信号同步移动）。

所以：问题是使用卡尔曼滤波器你的信号模型都是错误的。

更好的信号模型是在频域中：信号是低频的。噪音是高频的。因此，低通滤波器将在改进卡尔曼滤波器希望做的事情方面做得更好。

如果我将卡尔曼滤波器应用于您的问题（对您的方程进行一些修改以考虑和之间的差异），那么我得到下图。 $\mathbf{X}[k|k-1]$ $\mathbf{X}[k|k]$

带点的黑色曲线是measurement. 红色曲线是signal。蓝色曲线是卡尔曼滤波器的输出。绿色曲线是 5 点移动平均滤波器（简单的低通滤波器）的输出。

卡尔曼滤波器的误差平方和为118.63722。低通滤波器的相同数字是14.18046。显然，低通滤波器的信号模型拟合得更好，因为误差更小。

实现这一点的 R 代码如下。

# 26489
T <- 128;
t <- 0:(T-1)/T*70

signal <- 4.4 + sin(0.06*pi*t - 0.1) + sin(0.14*pi*t + 0.07) + sin(0.2*pi*t + 0.4)
noise <- sin(0.06 + 0.5*pi*t) + sin(0.1 - 0.68*pi*t) - sin(0.04 - 0.74*pi*t) + sin(0.03 - 0.93*pi*t)

measurement <- signal + noise


xkm1km1 <- matrix(c(4.4, 0),2,1)
Pkm1km1 <- matrix(c(1,0,0,1),2,2)
H <- matrix(c(1,0),1,2) 
A <- matrix(c(1,1,0,1),2,2)
Q <- 10^-6 * matrix(c(1.6,2.4,2.4,4.8),2,2)
R <- 10^-6 * matrix(c(2.35),1,1)

library("MASS") # For pseudo inverse ginv()

zhat <- t*0

for (k in 1:T)
{
  xkkm1 <- A %*% xkm1km1
  Pkkm1 <- A %*% Pkm1km1 %*% t(A) + Q
  K <- Pkkm1 %*% t(H) %*% ginv( H %*% Pkkm1 %*% t(H) + R)
  z <- matrix(c(measurement[k]), 1, 1)
  xkm1km1 <- xkkm1 + K %*% (z - H %*% xkkm1)
  Pkm1km1 <- (matrix(c(1,0,0,1),2,2) - K %*% H) %*% Pkkm1  
  zhat[k] <- as.numeric(H %*% xkkm1)
}

plot(t, measurement, type="b", pch=19)
lines(t, signal, col="red",  lwd=10)
lines(t, zhat, col="blue", lwd=5)

lpf <-  filter(measurement, c(0.2,0.2,0.2,0.2,0.2), circular = TRUE)
lines(t,lpf, col="green", lwd=4)

errs <- c(sum((signal - zhat)^2), sum((signal - lpf)^2))
print(errs)

其它你可能感兴趣的问题

上一篇导出 DTFTx [ 2 n ]x[2n] 下一篇滤波信号的方差