编辑：在考虑了更多之后，我想出了比锁相环简单得多的东西。

您遇到的问题是因为您正在使用棚车进行过滤。Boxcar 滤波器在频域中有很多波纹，因此如果您选择了错误的宽度，您将无法很好地衰减大约 10Hz 的信号。

如果您使用巴特沃斯滤波器，您将获得没有波纹的频率响应。 超过（比如说）5Hz 的所有东西都会衰减超过 -3db。巴特沃斯滤波器的计算成本也很低。

我去了http://www-users.cs.york.ac.uk/~fisher/mkfilter/trad.html并要求在 1000Hz 的采样率上使用截止频率为 5Hz 的三阶巴特沃斯低通滤波器并得到以下递归关系：

y[n] = (  1 * x[n- 3])
     + (  3 * x[n- 2])
     + (  3 * x[n- 1])
     + (  1 * x[n- 0])

     + (  0.9390989403 * y[n- 3])
     + ( -2.8762997235 * y[n- 2])
     + (  2.9371707284 * y[n- 1])

通过使用 Matlab 设计一个椭圆滤波器（可能是一个高于 3 阶的滤波器），您可能会得到更好的结果。您将在截止频率之外获得更锐利的衰减。

这是我关于锁相环的原始答案：

我会尝试离散时间锁相环。

所以是这样的：

这个想法是将您的输入信号乘以相同频率的正弦曲线，然后是低通滤波器。这会将输入信号的正弦部分转换为（几乎）零频率。通过观察低通滤波器输出的变化，您可以估计出您的估计频率与实际频率的差距。所以你反馈调整。

上图中乘法器的输出为：$$ Be^{2\pi i \theta_n} + Ae^{2 \pi i (\theta_n + \omega n)} + Ae^{2 \pi i (\ theta_n - \omega n)} + \varepsilon'(n)。$$ $\theta_n$ 是估计频率的运行总和，$\hat{\omega}_n$。当 $\hat{\omega}_n = \omega$ 时，$\theta_n-\omega n$ 是常数，所以 $A e^{2 \pi i (\theta_n - \omega n)}$ 是常数。

$$Be^{2\pi i \theta_n} + Ae^{2 \pi i (\theta_n + \omega n)} + Ae^{2 \pi i (\theta_n - \omega n)} + \varepsilon'(n).$$ $\theta_n$ is the running sum of the estimated frequency, $\hat{\omega}_n$. When $\hat{\omega}_n = \omega$ then $\theta_n-\omega n$ is constant, so $A e^{2 \pi i (\theta_n - \omega n)}$ is constant.

你想要的导数是角导数。给定来自低通滤波器的样本， $r_ne^{2\pi i \alpha_n}$ 和 $r_{n+1}e^{2\pi i \alpha_{n+1}}$ 你真的想要 $\alpha_ {n+1}-\alpha_{n}$。将低通样本 $n+1$ 除以低通样本 $n$ 并获取复杂的部分，您基本上得到了正确的结果。$r_ne^{2\pi i \alpha_n}$ and $r_{n+1}e^{2\pi i \alpha_{n+1}}$ you really want $\alpha_{n+1}-\alpha_{n}$. Dividing low-passed sample $n+1$ by low-passed sample $n$ and taking the complex part gets you essentially the right thing.

“调整 $\hat{\omega}_n$” 框也很重要。估计的相位误差可能有点嘈杂，因此您可能希望对相位误差进行低通滤波，并且您可能只想将 $\hat{\omega}_n$ 调整为估计相位误差的一小部分，而不是由整个相位误差。$\hat{\omega}_n$" box is also non-trivial. The estimated phase error is likely to be a little noisy so you might want to low-pass filter the phase error and you might want to only adjust $\hat{\omega}_n$ by a fraction of the estimated phase error, rather than by the entire phase error.

$\hat{\omega}_n$ is the estimated frequency, and $1/\hat{\omega}_n$ is the low-pass filter width you are looking for in your question. I think you can apply the same trick to the low-pass filter in the phase-locked loop. By dynamically adjusting the phase-locked loop low-pass filter width to $1/\hat{\omega}_n$ you completely eliminate $Be^{2\pi i \theta_n}$.

tl; dr：最后几段实际上给出了问题的可能答案，包括窗口化。

在您尝试改进方法的任何情况下，都值得看看有多少改进是可用的。这样你就可以决定改进的价值，决定它是否值得追求。

您试图避免的错误基本上是由信号中的正弦曲线的部分周期引入的。如果有完整数量的周期，那么取时间平均值是您将获得的最佳结果，因为正弦曲线在总和中消失了（我在这里假设零均值噪声）。

下图显示了信号中不完整周期将引入的误差量。 这显示了 1 秒内正弦信号（单位幅度）的平均值，其中 $f_s = 1$kHz。如果采样时间超过一秒钟，误差会减少（如您所见），误差将取决于正弦曲线的幅度。作为参考，经过两秒钟的采样，您最终会遇到这么多错误： $f_s = 1$kHz. If you sample for more than one second, the error will reduce (as you note), and the error will depend on the amplitude of the sinusoid. For reference, with two seconds of sampling, you end up with this much error:

您正确地注意到，您不一定能够使用 DFT 测量信号的确切频率。想到解决这个问题的一种方法可能是在从末端切掉不同数量的样本后计算信号的 FFT。这样做可以最大限度地减少频谱中的频谱泄漏（裙边）量，然后您将找到最接近正弦曲线完整周期数的近似值，并且可以假设信号均值大约是您将获得的最佳值.

您可以通过识别和去除频谱中的峰值，然后找到频谱中周围几个值的平均值来检测泄漏的最小化。当该平均值最小时（随着您继续从信号中删除样本，它再次开始上升），您就有了最佳位置。

编辑：忘了提一下，另一种简单的方法是将您的信号窗口化（例如，使用汉明窗口）。然后找到窗口信号的平均值，并将其放大以抵消窗口的影响。缩放值应该是 $$ \frac{N}{\sum_{n=1}^N w[n]} $$ 其中 $N$ 是信号中的样本数。尝试不同的窗口，看看什么效果最好。

$$\frac{N}{\sum_{n=1}^N w[n]}$$ where $N$ is the number of samples in your signal. Experiment with different windows to see what works best.

您的采样信号看起来像

$$x(n) = A\sin\left (\frac{2\pi fn}{f_s} + \phi\right ) + B + e(n)$$

其中$A$是幅度，$\phi$是相位，$f$是正弦信号的频率，$f_s$是采样频率，$B$是直流偏移，$e(n)$是一些加性噪声，希望是零均值。$A$、$f$、$\phi$、$B$ 和 $e(n)$ 未知，您想知道 $B$。$A$ is the amplitude, $\phi$ is the phase, $f$ is the frequency of the sinusoidal signal, $f_s$ is the sampling frequency, $B$ is the DC offset, and $e(n)$ is some additive noise, hopefully with zero mean. $A$, $f$, $\phi$, $B$, and $e(n)$ are unknown, and you would like to know $B$.

您对信号进行平均的想法是估计 $B$ 的一个完全有效的选择（如果加性噪声 $e(n)$ 的平均值确实为零）。您注意到的问题当然是您永远不知道采样数据中是否有整数个周期。要改进您的解决方案，您可以尝试以下方法。获取大量数据。你的块长度 $N$ 必须满足 $$N\gg \frac{f_s}{f}\tag{1}$$ 然后创建足够长度的子块（仍然满足（1））。要获得子块，只需从原始数据块中删除重叠块。通过这种方式，您可以相对于数据块的开头和结尾引入信号的随机相位，如果您对每个子块进行平均，然后对所有子块的结果进行平均，您将获得更好的结果总估计。$B$ (if the additive noise $e(n)$ has indeed a mean value of zero). The problem you have noted is of course that you never know if there is an integer number of periods in your sampled data. To improve on your solution you could try the following. Take a large block of data. Your block length $N$ must satisfy

$$N\gg \frac{f_s}{f}\tag{1}$$ Then create sub-blocks of sufficient length (still satisfying (1)). To get your sub-blocks simply cut out overlapping blocks from your original block of data. In this way you introduce a random phase of your signal with respect to the beginning and end of the data blocks, and if you average over each sub-block, and then average the results of all sub-blocks, you'll get a better total estimate. Say your estimate based on the $i^{th}$ sub-block is

$$\hat{B_{i}}=\frac{1}{M}\sum_{n=n_i}^{n_i+M-1}x(n)$$

$M$ 是子块中数据样本的数量，那么您对 ​​DC 偏移的最终估计将是$M$ the number of data samples in the sub-blocks, then your final estimate of the DC offset would be

$$\hat{B} = \frac{1}{L}\sum_{i=0}^{L-1}\hat{B_{i}}$$

其中 $L$ 是子块的数量。$L$ is the number of sub-blocks.