检测具有已知频率的正弦波的相位和幅度

Question

我从传感器收到正弦数据，其形式为(A + B(sin(n/N+a)))，其中 N 是样本总数，加上一些小噪声。我知道在N个样本（1000个样本）中，正弦波将完成一圈。信号如下所示：

我想使用尽可能少的数据来提取可变幅度“B”和相位“a”。换句话说，我想使用DSP尽快预测传感器的信号。我已经尝试过“相关性”，但没有用。

采用TMS320C000搭配FPU、TMU单元。

Answer 1

首先请注意，如果您的正弦波是周期性的N，则它实际上具有以下形式A+B*sin(2*pi*n/N + a)。

对于没有噪声且频率已知的信号，未知参数A、B和a可以通过最少 3 个样本获得。这可以通过首先求解以下线性方程组以获得B和来完成a：

然后使用回代得到A = x[0] - B*sin(a). 该解决方案可以通过以下代码实现：

double K = 2*PI/N;
// setup matrix
//   [sin(1/N)-sin(0/N) cos(1/N)-cos(0/N)]
//   [sin(2/N)-sin(1/N) cos(2/N)-cos(1/N)]
double a11 = sin(K);
double a12 = cos(K)-1;
double a21 = sin(2*K)-sin(K);
double a22 = cos(2*K)-cos(K);
// Compute 2x2 matrix inverse, and multiply by delta x vector
// see https://www.wolframalpha.com/input/?i=inverse+%7B%7Ba,+b%7D,+%7Bc,+d%7D%7D
double d   = 1.0/(a11*a22-a12*a21);
double y0  = d*(a22*(x[1]-x[0])-a12*(x[2]-x[1])); // B*cos(a)
double y1  = d*(a11*(x[2]-x[1])-a21*(x[1]-x[0])); // B*sin(a)
// Solve for parameters
double B   = sqrt(y0*y0+y1*y1);
double a   = atan2(y1, y0);
double A   = x[0] - B*sin(a);

由于您的信号包含一些噪声，因此使用更多样本对超定系统使用最小二乘近似解会获得更好的结果。这个超定系统可以写成：

具有以下定义：

解决方案如下：

然后，您需要在解决方案的准确性和所使用的样本数量之间进行权衡。对于使用示例的解决方案M，可以使用以下类似 C 的伪代码来实现：

// Setup initial conditions
double K   = 2*PI/N;
double s   = sin(K);
double c   = cos(K);
double sp  = s;
double cp  = c;
double sn  = s*cp + c*sp;
double cn  = c*cp - s*sp;
double t1 = s;
double t2 = c-1;
double b11 = 0.0;
double b12 = 0.0;
double b21 = 0.0;
double b22 = 0.0;
double z1  = 0.0;
double z2  = 0.0;
double dx  = 0.0;

// Iterative portion
for (int i = 0; i < M-1; i++)
{
  // B_{i,j} = (A^T A)_{i,j} = sum_{k} A_{k,i} A_{k,j}
  // For a 2x2 matrix B and a given "k", 
  // we have two values t1 & t2 which represent A_{k,1} and A_{k,2}
  b11 += t1*t1;
  b12 += t1*t2;
  b21 += t2*t1;
  b22 += t2*t2;

  // Update z_i = (A^T \Delta x)_i = \sum_k A_{k,i} (\Delta x)_i
  dx   = x[i+1]-x[i];
  z1  += t1 * dx;
  z2  += t2 * dx;

  // Update t1 & t2 for next term
  t1 = sn-sp;
  t2 = cn-cp;

  // Iteratively compute sin(2*pi*k/N) & cos(2*pi*k/N) using trig identities:
  //   sin(x+2pi/N) = sin(2pi/N)*cos(x) + cos(2pi/N)*sin(x)
  //   cos(x+2pi/N) = cos(2pi/N)*cos(x) - sin(2pi/N)*sin(x)
  sp = sn;
  cp = cn;
  sn = s*cp + c*sp;
  cn = c*cp - s*sp;
}

// Finalization

// Compute inverse of B
double dinv = 1.0/(b11*b22-b12*b21);
double binv11 =  b22*dinv;
double binv12 = -b21*dinv;
double binv21 = -b12*dinv;
double binv22 =  b11*dinv;

// Compute inv(B)*A^T \Delta x which gives us the 2D vector [B*cos(a) B*sin(a)]
double y1  = binv11*z1 + binv12*z2; // B*cos(a)
double y2  = binv21*z1 + binv22*z2; // B*sin(a)
// Solve for "B", "a" and "A" parameters
double B   = sqrt(y1*y1+y2*y2);
double a   = atan2(y2, y1);
double A   = x[0] - B*sin(a);

您可能会发现，在收到每个新样本时，对每个新样本执行一次循环迭代很方便。另外，如果您想在 和 上获得持续更新A，B并且a估计您只需要在每次迭代时执行终结部分（循环后的代码部分）。

最后，为了对输入峰值有一点额外的鲁棒性，您可以跳过b11、b12、b21、b22和的更新z1，z2以获取较大的dx：

dx = x[i+1]-x[i];
// either use fixed threshold (requires manual tuning) for simplicity
// or update threshold  dynamically as a fraction of B once a reasonable estimate of
// B is obtained. 
if (abs(dx) < threshold)
{
  b11 += t1*t1;
  b12 += t1*t2;
  b21 += t2*t1;
  b22 += t2*t2;

  z1  += t1 * dx;
  z2  += t2 * dx;
}
// always update t1, t2, sp, cp, sn, cn
...