我从传感器接收到一个正弦数据，格式为 (A + B(sin(n/N+a)))，其中 N 是样本总数，加上一些小噪音。我知道在 N 个样本（1000 个样本）中，正弦波将完成一圈。信号如下所示：
我想使用尽可能少的数据提取可变幅度“B”和相位“a”。换句话说，我想用DSP尽快预测传感器的信号。我已经尝试过“相关性”，但没有奏效。
将 TMS320C000 与 FPU、TMU 单元一起使用。
【问题讨论】：
首先请注意，如果您的正弦波是周期为N 的周期性正弦波，它实际上将是A+B*sin(2*pi*n/N + a) 的形式。
对于没有噪声和已知频率的信号，只需 3 个样本即可获得未知参数 A、B 和 a。这可以通过首先求解以下线性方程组得到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);
由于您的信号包含一些噪声，因此使用最小二乘法approximate solution to an over-determined system 使用更多样本可以获得更好的结果。这个超定系统可以写成：
具有以下定义：
然后给出解决方案：
然后，您需要在解决方案的准确性和使用的样本数量之间进行权衡。对于使用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
...
【讨论】：