Написал код для RLS фильтра для подавления шумов, но фильтр не работает. Где-то пол секунды идет трескотня, а затем тишина. Посмотрел в отладчике -переменные имеют значение 0.NAN, но где они могли переполнится? Привожу свой код. Прошу помощи
#define Ord (3) #define sigma (100000000) #define lambda (0.98f)
struct common_rls { float H[Ord]; float G[Ord]; float Ri[Ord*Ord]; // (Ri)tr = (Ri)* float X[Ord]; float Wrk[Ord]; }; struct common_rls rls;
void rls_init (void) { float diag=(float)100.0/sigma; int ind=0; volatile int i=0, j=0; for (i=0; i<Ord; i++) { rls.H[i] = 0.0; rls.G[i] =0.0; rls.X[i] =0.0; for (j=0; j<Ord; j++) { if (i==j) { rls.Ri[ind++] = diag; } else rls.Ri[ind++]=0; }
} } float rls_flt(float x) { // n - filter length // omega - forgetting factor // x - filter input // y - filter output // X - input delay line // G - Kalman gain, vector of length n // Ri - inverse of correlation matrix n x n, might be simmetrical packed // h - impulse response, vector of order n
volatile int i, j, ind=0; volatile float ssum, csum; volatile float error; volatile float y=0; static volatile float x_delay[2];
/*while (x!=0) {*/ for (i=1; i>0; i--) { x_delay[i] = x_delay[i-1]; } x_delay[0]=x; // error scalar // input delay line shift for (i=Ord-1; i>0; i--) { rls.X[i] = rls.X[i-1]; } rls.X[0] = x_delay[1]; error = x;//y; for (i=0; i<Ord; i++) { y=rls.H[i]*rls.X[i]; } error = x-y;
// adaptation // calculate Kalman Gain // ind=0 for (i=0; i<Ord; i++) { rls.Wrk[i] =0; for (j=0; j<Ord; j++) rls.Wrk[i] =rls.Wrk[i]+ rls.Ri[ind++]*(rls.X[j]); } csum = lambda; for (i=0; i<Ord; i++) { csum = csum + rls.X[i]*rls.Wrk[i]; } csum = 1/csum; for (i=0; i<Ord; i++) { rls.G[i] = rls.Wrk[i]*csum; } ssum = 1/lambda; ind =0; // calculate inverse correlation for (i=0; i<Ord; i++) {
for (j=0; j<Ord; j++) {
rls.Ri[ind] = (rls.Ri[ind] - rls.G[i]* (rls.Wrk[j]))*ssum; ind++; } } // filtering
// impulse response update for (i=0; i<Ord; i++) { rls.H[i] = rls.H[i] + rls.G[i]*error; }
/*}*/ return y; }
|