root/applications/pmsm/simulator_zdenek/ekf_example/mpf_double.cpp @ 1440

Revision 1440, 5.5 kB (checked in by smidl, 12 years ago)

rychly exp + vice clenu randn

RevLine 
[1380]1#include "mpf_double.h"
2#include "math.h"
3#include <stdio.h>
[1440]4#include "RArandn.h"
5#include "fast_exp.h"
[1380]6
[1439]7static floatx om[N];
8static floatx Pt[N];
9static floatx sth[N];
10static floatx cth[N];
11static floatx isdm[N];
12static floatx isqm[N];
13static floatx th[N];
14static floatx w[N];
15static floatx lwi[N];
[1380]16
[1439]17static floatx qth;
18static floatx qom;
19static floatx r;
[1380]20
[1439]21
22floatx qrandu() {
[1380]23        if ( randu_i==randu_last )
24                randu_i =&RArandu[0];
25        else
26                randu_i++;
27
28        return *randu_i;
29}
30
[1439]31floatx qrandn() {
[1380]32        if ( randn_i==randn_last )
33                randn_i =&RArandn[0];
34        else
35                randn_i++;
36
37        return *randn_i;
38}
39
40void resample() {
41        int N_babies[N];
[1439]42        floatx cumdist;
[1380]43        int i;
[1439]44        floatx ui;
45        int i_from=0;
46        int i_to=0;
[1380]47
[1439]48        floatx u0 = qrandu();
[1380]49
50        int j=0;
51        N_babies[0]=0;
52        cumdist = w[0];
53        for ( i = 0; i < N; i++ ) {
54                ui  = ( i + u0 ) / N;
55                while ( ui > cumdist ) {
56                        j++;
57                        N_babies[j]=0;
58                        cumdist+=w[j];
59                }
60                N_babies [j] ++;
61        }
62        while (j++<=(N-1)) { // delete all N_babies after j
63                N_babies[j]=0; 
64        }
65
66        // COPY
[1439]67
[1380]68        while ( i_from<N ) {
69                while ( N_babies[i_from]>1 ) { // 1 baby stays where it is
70                       
71                        while ( N_babies[i_to]>=1 ) i_to++; // find first empty slot
72
73                        //copy it
74                        Pt[i_to]=Pt[i_from];
75                        om[i_to]=om[i_from];
76                        th[i_to]=th[i_from];
77                        isdm[i_to]=isdm[i_from];
78                        isqm[i_to]=isqm[i_from];
79
80                        N_babies[i_to]++;
81                        N_babies[i_from]--;
82                }
83                i_from++;
84        }
85}
86
[1439]87void mpf_bayes ( const floatx isa, const floatx isb , const floatx usa, const floatx usb ) {
[1380]88
[1439]89        floatx isd;
90        floatx isq;
91        floatx usd;
92        floatx usq;
[1380]93
[1439]94        floatx Cd;
95        floatx Cq;
96        floatx CC;
97        floatx oCC;
[1380]98
[1439]99        floatx difid;
100        floatx difiq;
[1380]101
[1439]102        floatx zeta;
103        floatx Kd;
104        floatx Kq;
105        floatx ro;
106        floatx ypd;
107        floatx ypq;
108        floatx detRy;
109        floatx ydiffd;
110        floatx ydiffq;
111        floatx ydC;
[1380]112
[1439]113        floatx maxlwi, sumlwi;
[1380]114        int i;
[1439]115       
116        floatx *th_i, *Pt_i, *om_i, *isdm_i, *isqm_i;
117        floatx *w_i, *lw_i;
118        floatx *sth_i, *cth_i;
[1380]119// implementation starts here
[1439]120        th_i = &th[0];
121        Pt_i = &Pt[0];
122        om_i = &om[0];
123        isdm_i = &isdm[0];
124        isqm_i = &isqm[0];
125        sth_i = &sth[0];
126        cth_i = &cth[0];
[1380]127        for ( i=0; i<N; i++ ) {
128
[1439]129                // time update
130                *Pt_i = 1.0*1.0*(*Pt_i)+qom; // Pt is now predictive variance
131                *th_i += *om_i*_dt+qth *  qrandn() ;
132                while ( *th_i>M_PI ) *th_i -=2*M_PI;
133                while ( *th_i<-M_PI ) *th_i+=2*M_PI;
[1380]134
[1439]135                *sth_i=sin ( *th_i );
136                *cth_i=cos ( *th_i );
137               
138                isd = *cth_i*isa+*sth_i*isb;
139                isq = -*sth_i*isa+*cth_i*isb;
140                // process old voltage
141                usd = (*cth_i)*usa+(*sth_i)*usb;
142                usq = -(*sth_i)*usa+(*cth_i)*usb;
[1380]143
144                Cd = isq*_dt;
145                Cq = -_b - isd* _dt;
146
[1439]147                difid=isd- _ad  * (*isdm_i) - _cd *(usd);
148                difiq=isq- _aq  * (*isqm_i) - _cq *(usq);
[1380]149
150                CC=Cd*Cd+Cq*Cq;
[1439]151                zeta = *Pt_i/ ( r+*Pt_i*CC );
[1380]152                oCC = ( 1-zeta*CC );
153                ro = oCC/r;
154
[1439]155                Kd = *Pt_i*Cd*ro;
156                Kq = *Pt_i*Cq*ro;
[1380]157
[1440]158                (*Pt_i)*=( 1.0- ( Kd*Cd+Kq*Cq ) );
[1380]159
[1439]160                ypd = Cd**om_i;
161                ypq = Cq**om_i;
[1380]162
[1439]163                //detRy = ro/r;
164                *om_i +=  Kd* ( difid - ypd ) +Kq* ( difiq-ypq );
[1380]165                ydiffd = ( ypd-difid );
166                ydiffq = ( ypq-difiq );
167                ydC = ydiffd*Cd + ydiffq*Cq;
168
[1439]169                //lwi[i] = 0.5* ( log ( detRy ) + ( ydC*ydC*zeta -(ydiffd*ydiffd+ydiffq*ydiffq) ) /r ) ;
170                lwi[i] = 0.5*( ydC*ydC*zeta -(ydiffd*ydiffd+ydiffq*ydiffq) ) /r ;
[1380]171               
[1439]172                *isdm_i=isd;
173                *isqm_i=isq;
174               
175                // shift all counters
176                th_i++;
177                Pt_i++;
178                om_i++;
179                isdm_i++;
180                isqm_i++;
181                sth_i++;
182                cth_i++;
183
[1380]184        }
[1440]185/*      maxlwi=-1e10;
[1439]186        lw_i = &lwi[0];
187        for ( i=0;i<N;i++,lw_i++ ) {
188                if ( *lw_i>maxlwi ) maxlwi=*lw_i;
[1440]189        }*/
[1439]190        lw_i = &lwi[0];
191        w_i = &w[0];
192        for ( i=0;i<N;i++,lw_i++,w_i++ ) {
[1440]193//              *lw_i-=maxlwi;
194                (*w_i)=EXP ( *lw_i ); // always resample -> wmin=1/n;
[1439]195        }
[1380]196
197        sumlwi=0.0;
[1439]198        w_i = &w[0];
199        for ( i=0;i<N;i++,w_i++ ) sumlwi+=*w_i;
200        w_i = &w[0];
201        for ( i=0;i<N;i++,w_i++ ) *w_i/=sumlwi;
202       
203        // NEFF
204        floatx sw2=0.0;
205        for ( i=0;i<N;i++ ) sw2+=w[i]*w[i];
206        printf("%f\n",1./sw2*N);
207       
[1380]208        resample();
209}
210
211
[1439]212void mpf_init(floatx qom0, floatx qth0, floatx r0) {
213        int i;
[1380]214        r=r0;
215        qth=qth0;
216        qom=qom0;
[1439]217
[1380]218        for ( i=0; i<N; i++ ) {
[1439]219                th[i]=qrandu() *2*M_PI - M_PI;
220                sth[i]=sin ( th[i] );
221                cth[i]=cos ( th[i] );
222               
[1380]223                om[i]=0;
[1439]224                Pt[i]=10.;
[1380]225                isdm[i]=0;
226                isqm[i]=0;
227                w[i]=1.0/N;
228        }
229
230}
[1439]231void mpf_mean ( floatx *Ecosth, floatx *Esinth, floatx *Eome ) {
[1380]232        int i;
[1440]233        floatx *w_i=&w[0];
234        floatx *cth_i=&cth[0];
235        floatx *sth_i=&sth[0];
236        floatx *om_i=&om[0];
[1381]237        *Ecosth=0.0;
238        *Esinth=0.0;
239        *Eome=0.0;
[1440]240        for ( i=0;i<N;i++,w_i++,cth_i++,sth_i++,om_i++ ) {
241                *Ecosth+=(*w_i)*(*cth_i);
242                *Esinth+=(*w_i)*(*sth_i);
243                *Eome+=(*w_i)*(*om_i);
[1380]244        }
245
246}
[1439]247void mpf_th(floatx th1[N]){
248        int i;
249        for (i=0;i<N;i++) th1[i]=th[i];
[1380]250}
[1438]251
[1439]252floatx kalman_om( const floatx isa, const floatx isb , const floatx usa, const floatx usb, const floatx th ) {
[1438]253
[1439]254        floatx isd;
255        floatx isq;
256        floatx usd;
257        floatx usq;
[1438]258
[1439]259        floatx Cd;
260        floatx Cq;
261        floatx CC;
262        floatx oCC;
[1438]263
[1439]264        floatx difid;
265        floatx difiq;
[1438]266
[1439]267        floatx zeta;
268        floatx Kd;
269        floatx Kq;
270        floatx ro;
271        floatx ypd;
272        floatx ypq;
273        floatx detRy;
274        floatx ydiffd;
275        floatx ydiffq;
276        floatx ydC;
[1438]277
[1439]278        static floatx Pt;
279        static floatx om;
280        floatx cth, sth;
281        static floatx isdm, isqm;
[1438]282
283       
284                Pt = 1.0*1.0*Pt+qom; // Pt is now predictive variance
285                //while ( th>M_PI ) th=th-2*M_PI;
286                //while ( th<-M_PI ) th=th+2*M_PI;
287
288                sth=sin ( th );
289                cth=cos ( th );
290
291                isd = cth*isa+sth*isb;
292                isq = -sth*isa+cth*isb;
293                usd = cth*usa+sth*usb;
294                usq = -sth*usa+cth*usb;
295
296                Cd = isq*_dt;
297                Cq = -_b - isd* _dt;
298
299                difid=isd- _ad  *isdm - _cd *usd;
300                difiq=isq- _aq  *isqm - _cq *usq;
301
302                CC=Cd*Cd+Cq*Cq;
303                zeta = Pt/ ( r+Pt*CC );
304                oCC = ( 1-zeta*CC );
305                ro = oCC/r;
306
307                Kd = Pt*Cd*ro;
308                Kq = Pt*Cq*ro;
309
310                Pt=Pt* ( 1- ( Kd*Cd+Kq*Cq ) );
311
312                ypd = Cd*om;
313                ypq = Cq*om;
314
315                detRy = ro/r;
316                om = om + Kd* ( difid - ypd ) +Kq* ( difiq-ypq );
317                ydiffd = ( ypd-difid );
318                ydiffq = ( ypq-difiq );
319                ydC = ydiffd*Cd + ydiffq*Cq;
320
321                isdm=isd;
322                isqm=isq;
323               
324                return om;
325}
Note: See TracBrowser for help on using the browser.