信息处理 - 估计具有部分已知频率的谐波信号的线性组合参数 - 吾爱随笔录

估计具有部分已知频率的谐波信号的线性组合参数

信息处理离散信号窗函数参数估计

2022-02-02 00:25:28

考虑以下形式的真实信号：

\begin{aligned} x (n) = & \sum_{m = 1}^{M} (c_{m} \cos (2 π f_{m} n) + s_{m} \sin (2 π f_{m} n)) = \\ \sum_{m = 1}^{M} a_{m} \cos (2 π f_{m} n + ϕ_{m}) \end{aligned}

$\begin{align*} x(n) = &\sum_{m=1}^{M} \left( c_m \cos(2\pi f_m n) + s_m \sin(2\pi f_m n) \right) =\\ &\sum_{m=1}^{M} a_m \cos(2\pi f_m n + \phi_m) \end{align*}$ 其中的信号值

x (n)

$x(n)$ 以

n = 1, \dots, N >> M

$n = 1,\dots, N >> M$ 和一些频率

f_{m}

$f_m$ 是已知的实数

(0.0, 0.5)

$(0.0,0.5)$ .

最好的估算方法是什么 $c_m$ 和 $s_m$ 实际参数？

估计形式：

c_{m} + j s_{m} = 2 \frac{\sum_{n = 1}^{N} x (n) \exp (2 π j f_{m} n)}{N}

$c_m + j s_m = 2 \frac{\displaystyle\sum_{n=1}^{N} x(n) \exp(2 \pi j f_m n)}{N}$ 不是很准确（

j^{2} = - 1

$j^2=-1$ )，但加窗（例如余弦窗）显着提高了估计精度：

c_{m} + j s_{m} = 2 \frac{\sum_{n = 1}^{N} w (n) x (n) \exp (2 π j f_{m} n)}{\sum_{n = 1}^{N} w (n)}

$c_m + j s_m = 2 \frac{\displaystyle\sum_{n=1}^{N} w(n) x(n) \exp(2 \pi j f_m n)}{\displaystyle\sum_{n=1}^{N} w(n)}$

为什么窗口在这里如此有效？它是某种众所周知的公式吗？

如果所有项的频率都是已知的，最小二乘法也能得到很好的结果。

在并非所有频率都已知的情况下，还有其他好的方法吗？也许基于 SVD 或某种正交化。

Mathematica 中的测试代码：

F1 = 0.123000 ;
F2 = 0.312874 ;

LENGTH = 1024 ;
RANGE = Range[LENGTH] ;

SIGNAL = 0.0 ;
SIGNAL = SIGNAL + 1.0*Cos[2*Pi*F1*RANGE] + 0.5*Sin[2*Pi*F1*RANGE] ;
SIGNAL = SIGNAL + 0.05*Cos[2*Pi*F2*RANGE] + 0.01*Sin[2*Pi*F2*RANGE] ;

WINDOW = HannWindow[N[Rescale[Range[LENGTH],{1,LENGTH},{-1,1}]]]^2 ;

(* NO WINDOW *)
2*(Total[SIGNAL*Exp[2*Pi*I*F1*RANGE]]/LENGTH)               
2*(Total[SIGNAL*Exp[2*Pi*I*F2*RANGE]]/LENGTH)               
(* OUTPUT: 0.9994929536863852`\[VeryThinSpace]+0.49994639046951617` I *)
(* OUTPUT: 0.04819835064062154`\[VeryThinSpace]+0.012127390241918076` I *)

(* WITH WINDOW *)
2*(Total[SIGNAL*WINDOW*Exp[2*Pi*I*F1*RANGE]]/Total[WINDOW]) 
2*(Total[SIGNAL*WINDOW*Exp[2*Pi*I*F2*RANGE]]/Total[WINDOW])
(* OUTPUT: 0.9999999999775656`\[VeryThinSpace]+0.5000000000025541` I *)
(* OUTPUT: 0.05000000004559643`\[VeryThinSpace]+0.009999999958900167` I *)

(* LEAST SQUARES *)
(* ONE KNOWN *)
LeastSquares[{Cos[2*Pi*F1*RANGE],Sin[2*Pi*F1*RANGE]}//Transpose,SIGNAL]
(* ALL KNOWN *)
LeastSquares[{Cos[2*Pi*F1*RANGE],Sin[2*Pi*F1*RANGE],Cos[2*Pi*F2*RANGE],Sin[2*Pi*F2*RANGE]}//Transpose,SIGNAL]
(* OUTPUT: {0.999957445314386`,0.4999499273288035`} *)
(* OUTPUT: {1.`,0.49999999999999983`,0.05000000000000022`,0.010000000000000056`} *)

2个回答

底线：OP 正在经历频谱泄漏的影响以及使用离散傅里叶变换进行窗口化的主要动机。加窗方法获得更好的估计精度的原因解释如下：

OP 的方程是真实正弦曲线的总和，他确定这些正弦曲线的估计器是离散傅里叶变换 (DFT) 的一种形式。如果没有额外的窗口（这是一个矩形窗口），众所周知，DFT 会遭受大量的频谱泄漏，因为矩形窗口的内核（离散时间傅里叶变换或矩形窗口的 DTFT）混叠 Sinc 函数的形式）在频率上与每个真实正弦波的基本两个脉冲频率进行卷积。由于混叠 Sinc 函数随着频率偏离代表音调的脉冲位置而相对逐渐减小，并且具有如此高的旁瓣，许多其他频率区间受到影响，导致它们自己的测量出现错误（以及动态范围显着降低或在存在强音时检测较低电平音调的能力）。选择的窗口具有低得多的旁瓣电平，并且频率滚降快得多，因此任何给定的音调对存在的其他音调的估计影响要小得多。

如果我们观察 DTFT 而不是 DFT，这一点会更清楚，因为 DTFT 是频率的连续函数，而 DFT 是该连续函数上的样本。通过观察 DTFT，我们可以更轻松地全面直观地了解上述内容。此外，我们还看到了如何逼近 DTFT（通过在进行 DFT 之前进行零填充）也可以用来逼近频率）。我用 OP 的原始频率图展示了这一点，我们可以直接看到窗口的效果：

首先考虑 OP 的 1024 个样本示例，仅 F1 单独使用和不使用额外的窗口。首先没有额外的窗口（矩形窗口） DTFT 的对数幅度图与对应的 F1 的 1024 个 DFT 样本如下图所示：

这是专门的 DTFT

s i g 1 = 1.0 \cos (2 π f_{1} n) + 0.5 \sin (2 π f_{1} n)

$sig1 = 1.0\cos(2\pi f_1 n) + 0.5\sin(2\pi f_1 n)$

对于 n = 0 到 1023。而 f1 = 0.123

（注：以上只是单音 $f1$ 具有额外的相移，可以从关系中得出 $a\cos(x)+b\cos(x) = \sqrt{a^2+b^2}\cos(x-\tan^{-1}(b/a))$

仅 F1 信号的幅度谱 - 矩形窗口

如果我们放大其中一个光谱峰值，我们可以清楚地看到 Sinc 形内核（在这种情况下，考虑到样本数量，混叠将非常小，因此这是实际 Sinc 的非常接近的近似值）与卷积F1 频率位置的潜在脉冲，我们可以看到频谱泄漏以及它如何影响其他频率位置。F1 为 0.123 非常接近 bin 编号 126，如下所示 $1024 \times .123 = 125.952$ .

放大 F1 峰值位置

下面我展示了 OP 使用的两个频率分量，但仍然作为单独的信号来演示每个的频谱泄漏如何影响另一个的测量。

这是具体的

s i g 1 = 1.0 \cos (2 π f_{1} n) + 0.5 \sin (2 π f_{1} n)

$sig1 = 1.0\cos(2\pi f_1 n) + 0.5\sin(2\pi f_1 n)$

s i g 2 = 0.05 \cos (2 π f_{2} n) + 0.01 \sin (2 π f_{2} n)

$sig2 = 0.05\cos(2\pi f_2 n) + 0.01\sin(2\pi f_2 n)$

对于 n = 0 到 1023 和 f1 = 0.123 和 f2 = 0.312874

F1 和 F2 信号的幅度谱（单独） - 矩形窗口

请注意，如果我们放大 F2 的峰值，来自 F1 的泄漏有多显着将影响 F2 幅度的估计，因为 F1 的泄漏仅低 29 dB！如果 F2 的幅度降低，则产生的误差会更加明显。（另请注意，我们如何使用 DTFT 将 F2 频率精确定位在其分数 bin 位置 $1024 \times 0.312874 \approx 320.3830$ 这是下图中橙色曲线的峰值所在的位置，但一旦与蓝色曲线给出的 F2 信号结合，这也会导致频率误差，因为蓝色曲线的峰值稍微偏右橙色 F1 曲线。

放大 F2 峰位置 - 使用矩形窗口

如果我们然后计算具有 F1 和 F2 的复合波形的 DTFT，我们会看到在 F2 峰值位置附近产生的幅度和频率误差，如下图所示（其中绿色是复合波形的 DTFT 和 F1和 F2 组件）。

该误差由两个波形的相干和预测，因此可以从矩形窗口内核和两个信号的相对幅度和频率偏移的知识中预测。鉴于波形很复杂，分量的相位会影响求和结果（因此干扰可能会降低或增加受影响分量的幅度）。

在此示例中，以下是 F2 峰值位置的信号贡献（不是上面的 dB 标度以直接显示总和）。

F2 单独在 bin 320.38：幅度 = 52.10（给出为 52.00 -j3.26）

F1 单独在 bin 320.38：幅度 = 1.70（给出为 1.42-j0.94）

F1 + F2 在 bin 320.38：幅度 = 53.58（给出为 53.42-j4.19）

其中，测量的 F1 + F2 仅通过添加复杂分量时的 F1 和 F2 之和来预测。

最后，如果我们使用加窗观察上面的图（我们将使用与 OP 相同的 Hann 窗），我们可以看到频谱泄漏的显着减少，以及这将如何有助于更好地估计频率和基础信号的幅度。

F1 和 F2 信号的幅度谱（单独）- Hann Window

放大 F2 频谱峰值现在显示 F1 的泄漏现在降低了 100 dB 以上，但我们也看到了窗口化的缺点 - 内核的主瓣更宽导致频率选择性降低。（矩形窗口的选择性最好，主瓣最窄，但由于旁瓣电平高且滚降慢，动态范围最差）。

放大 F2 峰位置 - 使用 Hann 窗口

有关 DTFT 与 DFT 以及使用零填充近似 DTFT 的更多信息，请参阅从 DFT 样本中获取 DTFT和对于 2D 信号，可以说频率响应与傅里叶变换相同吗？

有关窗口化的更多详细信息，请参阅：为什么要使用 Hann 或 Bartlett 窗口？

有关频谱泄漏的更多详细信息，请参阅：FFT 中旁瓣的直觉

高效加工的建议

假设目标是只找到每个音调的幅度，如果 0.01 dB 的精度就足够了，建议使用平顶窗口，它提供近 100 dB 的旁瓣抑制，同时在向外延伸的频率中保持非常平坦的主瓣到 $\pm 0.5$ 垃圾箱。这意味着不需要进一步的插值来估计位于分数频段频率的音调的幅度，以获得在实际幅度的 0.01 dB 范围内的估计值。

下面显示了一个只有 128 个样本的平顶窗口的内核 (DTFT) 示例，显示在 5 个 bin 偏移处达到了近 100 dB 的旁瓣抑制。

放大以查看平坦度，我们看到在 128 个点上，我们保持小于 0.01 dB 的平坦度 $\pm 0.5$ 垃圾箱。

如果 $0.01$ dB 幅度精度是可以接受的，最小化处理的解决方案是选择所需的最大 DFT 大小，以确保基于最小音调接近度没有两个音调间隔小于 5 个 bin（在大约 4.6 个 bin 处发生的 58.7 dB 泄漏会导致最大幅度误差为 0.01 dB）。这将设置要处理的最小样本数，之后可以乘以平顶窗口并在最大箱足以提供幅度估计的地方进行fft。

这是一篇有趣的论文，其中包含有关根据任何通带误差规范（以牺牲旁瓣性能为代价）设计平顶窗口以及其他优化平顶解决方案的更多信息：使用半无限优化的平顶窗口的设计和改进。

好吧，这些结果比我预期的差一点。我要再看一遍。

这是信号发生器：

将 numpy 导入为 np

#=================================================== ==============================
定义主（）：

        N = 64
        
        x = np.zeros(N)

        对于范围内的n（N）：
          x[n] = 1.0 * np.cos( 2 * np.pi * 0.123 * n ) \
               + 0.5 * np.sin（2 * np.pi * 0.123 * n）\
               + 0.1 * np.cos（2 * np.pi * 0.378 * n）\
               + 0.0001 * np.sin(2 * np.pi * 0.378 * n)

        f = open("TwoTone.txt", "w")
        
        对于范围内的n（N）：
          f.write(str(x[n])+“\n”)
                       
        f.close()               
                  
#=================================================== ==============================
主要的（）

这是16点案例：

16
峰值为 2 657.100
峰值 6 65.527

    2 1.967574268 1.127686414 -0.462437761 9.677939450
    6 6.096984210 0.112391665 -0.215513086 8.704864641
theShiftSquareSum 1.18480464437

    2 1.967574268 1.118934768 -0.461595191 0.000064761
    6 6.096984210 0.113199253 -0.209107938 0.047386657
theShiftSquareSum 0.000153884710745

    2 1.967574268 1.118950503 -0.461589923 0.000002945
    6 6.096984210 0.113177001 -0.209377985 0.047165809
theShiftSquareSum 1.11779322854e-07

    2 1.967574268 1.118950601 -0.461589890 0.000002946
    6 6.096984210 0.113177942 -0.209366691 0.047171378
theShiftSquareSum 1.967368058e-10

    2 1.967574268 1.118950602 -0.461589890 0.000002946
    6 6.096984210 0.113177902 -0.209367163 0.047171141
theShiftSquareSum 3.44244306479e-13

    2 1.967574268 1.118950602 -0.461589890 0.000002946
    6 6.096984210 0.113177904 -0.209367144 0.047171151
theShiftSquareSum 6.02574348723e-16

    2 1.968 1.119 -0.462 0.12297 1.00185 0.49835
    6 6.097 0.113 -0.209 0.38106 0.11071 0.02352

这是 64 点的案例：

64
高峰在 8 1596.112
峰值为 24 152.103

    8 7.871922908 1.125979444 -0.458494786 215.881795607
   24 24.200735389 0.108911265 -0.137066216 196.820341835
移位平方和 1.76424051554

    8 7.871922908 1.118802923 -0.461915642 0.000058885
   24 24.200735389 0.109116108 -0.137631501 0.048358225
theShiftSquareSum 5.61078445351e-05

    8 7.871922908 1.118765771 -0.461961717 0.000005558
   24 24.200735389 0.109117039 -0.137631515 0.048362634
theShiftSquareSum 3.29842827826e-09

    8 7.871922908 1.118765273 -0.461962339 0.000005549
   24 24.200735389 0.109117039 -0.137631515 0.048362634
theShiftSquareSum 5.85271671495e-13

    8 7.871922908 1.118765266 -0.461962348 0.000005549
   24 24.200735389 0.109117039 -0.137631515 0.048362634
theShiftSquareSum 1.06705350002e-16

    8 7.872 1.119 -0.462 0.12300 1.00150 0.49864
   24 24.201 0.109 -0.138 0.37814 0.10809 0.01497

最后三列应与参数匹配。我期望的准确性至少 $10^{-15}$ 在这种无噪音的情况下。

其它你可能感兴趣的问题

上一篇如何区分 FFT 功率泄漏或关闭频率下一篇在某些采样率下采样低频声音不起作用