Erkennung von Phase und Amplitude eines Sinus mit bekannter Frequenz

Question

Ich erhalte sinusförmige Daten von einem Sensor in der Form (A + B (sin (n/N + a))), wobei N die Summe ist Anzahl der Samples plus ein paar kleine Geräusche. Ich weiß, dass bei N Samples (1000 Samples) die Sinuswelle eine Umdrehung vollendet. Das Signal sieht wie folgt aus:

ich die Variable Amplitude "B" und Phase "a", mit so wenig Daten wie möglich extrahieren möchten. Mit anderen Worten, ich möchte das Signal des Sensors so bald wie möglich mit DSP vorhersagen. Ich habe bereits "Korrelation" versucht, aber es hat nicht funktioniert.

Verwenden von TMS320C000 mit FPU, TMU-Einheit.

Answer 1

Beachten Sie zuerst, dass wenn Ihre Sinuswelle periodisch mit der Periode N ist, sie tatsächlich die Form A+B*sin(2*pi*n/N + a) hätte. Für ein Signal ohne Rauschen und einer bekannten Frequenz konnten die unbekannten Parameter A, B und a mit nur 3 Proben erhalten werden. Dies kann durch zunächst Lösen der folgenden linearen Gleichungssystems durchgeführt werden, um zu erhalten und Ba: Substitution

dann unter Verwendung von Rück A = x[0] - B*sin(a) zu erhalten. Diese Lösung kann mit dem folgenden Code implementiert werden:

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); 

Da Ihr Signal enthält einige Geräusche, würden Sie bessere Ergebnisse mit einem der kleinsten Quadrate approximate solution to an over-determined system die Nutzung von mehr Proben erhalten.

Mit den folgenden Definitionen:

die Lösung wird dann gegeben als:

du überbestimmtes System kann geschrieben werden Dann muss ein Kompromiss zwischen der Genauigkeit der Lösung und der Anzahl gefunden werden der Probe verwendet. Für eine Lösung M Proben verwendet wird, kann dies die folgenden C-ähnlichen Pseudocode implementiert werden:

// 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); 

Sie finden es praktisch eine Wiederholung der Schleife für jede neue Probe auszuführen, wie Sie sie empfangen. Wenn Sie fortlaufende Updates für Ihre A, B und a Schätzung erhalten möchten, müssen Sie lediglich den Finalisierungsteil (den Teil des Codes nach der Schleife) bei jeder Iteration ausführen.

schließlich für ein wenig zusätzliche Robustheit Eingang Spikes, können Sie Updates auf b11, überspringen b12, b21, b22, z1 und z2 für große 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 
...