真正的零均值高斯白噪声，独立于一个干净的信号$x$和已知的方差被添加到$x$产生一个有噪声的信号$ y。：
$$Y_k = \frac{1}{N}\sum_{n=0}^{N-1}e^{-i2\pi kn/N}y_n.\tag{1}$$
这仅用于上下文，我们将在频域中定义噪声方差，因此归一化（或缺乏归一化）并不重要。时域的高斯白噪声就是频域的高斯白噪声，参见问题：“什么是高斯白噪声的离散傅里叶变换的统计量？ ”。因此我们可以写：
$$Y_k = X_k + Z_k,$$
其中$X$和$Z$是干净信号和噪声的 DFT，$Z_k$是遵循圆对称复高斯方差分布$\sigma^2$的噪声箱。$Z_k$的实部和虚部各自独立地遵循方差$\frac{1}{2}\sigma^2$的高斯分布。我们将 bin $Y_k$的信噪比 (SNR)定义为：
$$\mathrm{SNR} = \frac{\sigma^2}{|X_k|^2}.$$
然后尝试通过频谱减法来降低噪声，从而在保持原始相位的同时独立地减小每个 bin $Y_k$的幅度（除非在幅度减小中 bin 值变为零）。减少形成了对干净信号的 DFT 的每个 bin 的绝对值的平方$|X_k|^2$的估计$\widehat{|X_k|^2}$ ：
$$\widehat{|X_k|^2} = |Y_k|^2 - \sigma^2,\tag{2}$$
其中$\sigma^2$是每个 DFT bin 中已知的噪声方差。为简单起见，我们不考虑甚至$N$的$k = 0,$或$k = N/2$，这是真实$x.$的特殊情况。在低 SNR 下，(2) 中的公式有时可能结果为负$\widehat{| X_k|^2}.$我们可以通过从下面将估计值限制为零来消除这个问题，重新定义：
$$\widehat{|X_k|^2} = \max\left(|Y_k|^2 - \sigma^2,\,0\right).\tag{3}$$
图 1. 样本大小为$10^5,$的 Monte Carlo 估计： 实心：通过$\widehat{|X_k|}$估计$|X_k|$与使用$估计相比，误差平方和的增益|Y_k|,$虚线：通过$\widehat{|X_k|^2}$估计 $|X_k|^ 2 $ 与使用$|Y_k|^2,$进行估计相比
，平方误差的增益：用$\widehat{|X_k|}e^{i\arg(Y_k)}$估计 $X_k$ 与用$ Y_k估计相比，平方误差和的增益。 $\widehat{ | X_k的定义使用来自 (3) 的|^2}$ 。
问题：是否有另一个$|X_k|$或$|X_k|^2$的估计值在不依赖$Y_k$分布的情况下改进了（2）和（3） ？
我认为这个问题相当于用已知参数$\sigma_\mathrm{Rice} = \frac{估计Rice 分布（图 2）的参数$\displaystyle{\nu_\mathrm{Rice}}$的平方\sqrt{2}}{2}\sigma,$给定一个观察值。
图 2. Rice 分布是从一个点到原点的距离$R$的分布，该点遵循二元循环对称正态分布，均值绝对值为$\nu_\mathrm{Rice}，$方差$2\sigma_\mathrm {Rice}^2 = \sigma^2$和分量方差$\sigma_\mathrm{Rice}^2 = \frac{1}{2}\sigma^2.$
我发现了一些似乎相关的文献：
- Jan Sijbers、Arnold J. den Dekker、Paul Scheunders 和 Dirk Van Dyck，“Rician 分布参数的最大似然估计”，IEEE Transactions on Medical Imaging（第 17 卷，第 3 期，1998 年 6 月）（doi，pdf）。
用于估计曲线的 Python 脚本 A
可以扩展此脚本以在答案中绘制估计曲线。
import numpy as np
from mpmath import mp
import matplotlib.pyplot as plt
def plot_est(ms, est_as):
fig = plt.figure(figsize=(4,4))
ax = fig.add_subplot(1, 1, 1)
if len(np.shape(est_as)) == 2:
for i in range(np.shape(est_as)[0]):
plt.plot(ms, est_as[i])
else:
plt.plot(ms, est_as)
plt.axis([ms[0], ms[-1], ms[0], ms[-1]])
if ms[-1]-ms[0] < 5:
ax.set_xticks(np.arange(np.int(ms[0]), np.int(ms[-1]) + 1, 1))
ax.set_yticks(np.arange(np.int(ms[0]), np.int(ms[-1]) + 1, 1))
plt.grid(True)
plt.xlabel('$m$')
h = plt.ylabel('$\hat a$')
h.set_rotation(0)
plt.show()
图 1 的 Python 脚本 B
此脚本可以扩展为答案中的误差增益曲线。
import math
import numpy as np
import matplotlib.pyplot as plt
def est_a_sub_fast(m):
if m > 1:
return np.sqrt(m*m - 1)
else:
return 0
def est_gain_SSE_a(est_a, a, N):
SSE = 0
SSE_ref = 0
for k in range(N): #Noise std. dev = 1, |X_k| = a
m = abs(complex(np.random.normal(a, np.sqrt(2)/2), np.random.normal(0, np.sqrt(2)/2)))
SSE += (a - est_a(m))**2
SSE_ref += (a - m)**2
return SSE/SSE_ref
def est_gain_SSE_a2(est_a, a, N):
SSE = 0
SSE_ref = 0
for k in range(N): #Noise std. dev = 1, |X_k| = a
m = abs(complex(np.random.normal(a, np.sqrt(2)/2), np.random.normal(0, np.sqrt(2)/2)))
SSE += (a**2 - est_a(m)**2)**2
SSE_ref += (a**2 - m**2)**2
return SSE/SSE_ref
def est_gain_SSE_complex(est_a, a, N):
SSE = 0
SSE_ref = 0
for k in range(N): #Noise std. dev = 1, X_k = a
Y = complex(np.random.normal(a, np.sqrt(2)/2), np.random.normal(0, np.sqrt(2)/2))
SSE += abs(a - est_a(abs(Y))*Y/abs(Y))**2
SSE_ref += abs(a - Y)**2
return SSE/SSE_ref
def plot_gains_SSE(as_dB, gains_SSE_a, gains_SSE_a2, gains_SSE_complex, color_number = 0):
colors = plt.rcParams['axes.prop_cycle'].by_key()['color']
fig = plt.figure(figsize=(7,4))
ax = fig.add_subplot(1, 1, 1)
if len(np.shape(gains_SSE_a)) == 2:
for i in range(np.shape(gains_SSE_a)[0]):
plt.plot(as_dB, gains_SSE_a[i], color=colors[i], )
plt.plot(as_dB, gains_SSE_a2[i], color=colors[i], linestyle='--')
plt.plot(as_dB, gains_SSE_complex[i], color=colors[i], linestyle=':')
else:
plt.plot(as_dB, gains_SSE_a, color=colors[color_number])
plt.plot(as_dB, gains_SSE_a2, color=colors[color_number], linestyle='--')
plt.plot(as_dB, gains_SSE_complex, color=colors[color_number], linestyle=':')
plt.grid(True)
plt.axis([as_dB[0], as_dB[-1], 0, 2])
plt.xlabel('SNR (dB)')
plt.ylabel('SSE gain')
plt.show()
as_dB = range(-40, 41)
as_ = [10**(a_dB/20) for a_dB in as_dB]
gains_SSE_a_sub = [est_gain_SSE_a(est_a_sub_fast, a, 10**5) for a in as_]
gains_SSE_a2_sub = [est_gain_SSE_a2(est_a_sub_fast, a, 10**5) for a in as_]
gains_SSE_complex_sub = [est_gain_SSE_complex(est_a_sub_fast, a, 10**5) for a in as_]
plot_gains_SSE(as_dB, gains_SSE_a_sub, gains_SSE_a2_sub, gains_SSE_complex_sub, 1)
