既知の周波数を有するサインの位相と振幅

Question

（A + B（sin（n/N + a）））の形のセンサから正弦波のデータを受け取ります。ここでNは合計サンプル数、いくつかの小さなノイズ。私はN個のサンプル（1000サンプル）で、正弦波が1回転することを知っています。私は、できるだけデータを使用して、可変振幅「B」と位相「A」を抽出したい

：信号はこのようになります。換言すれば、DSPを使用してできるだけ早くセンサの信号を予測したいと考えています。私はすでに "相関"を試みたが、うまくいかなかった。

FPU、TMUユニットでTMS320C000を使用する。

Answer 1

まず、正弦波が期間Nの周期的である場合、実際にはA+B*sin(2*pi*n/N + a)の形式になることに注意してください。

ノイズがなく、既知の周波数の信号の場合、未知のパラメータA,Bおよびaは、わずか3サンプルで取得できました。

その後、A = x[0] - B*sin(a)を得るために、置換を使用してバック：これは最初のBと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のような擬似コードを使用して実装することができます：あなたは、あなたがそれらを受け取るようが便利それぞれの新しいサンプルに対してループの1回の繰り返しを実行するために見つけることができ

// 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の更新を継続したい場合は、繰り返しごとにファイナライゼーション部分（ループの後のコード部分）を実行するだけで済みます。

最後に、入力スパイクに余分な堅牢性の少しのために、あなたは大規模なdxためb11、b12、b21、b22、z1とz2の更新をスキップすることができます：

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 
...