hige dura
Dependencies: FatFileSystem mbed
main.cpp@0:813b917360ac, 2012-04-02 (annotated)
- Committer:
- higedura
- Date:
- Mon Apr 02 12:22:30 2012 +0000
- Revision:
- 0:813b917360ac
- Child:
- 1:8048a8bcde59
Who changed what in which revision?
| User | Revision | Line number | New contents of line |
|---|---|---|---|
| higedura | 0:813b917360ac | 1 | // dot-HR EKF |
| higedura | 0:813b917360ac | 2 | |
| higedura | 0:813b917360ac | 3 | #include "mbed.h" |
| higedura | 0:813b917360ac | 4 | #include "SDFileSystem.h" |
| higedura | 0:813b917360ac | 5 | #include "ADXL345_I2C.h" |
| higedura | 0:813b917360ac | 6 | #include "ITG3200.h" |
| higedura | 0:813b917360ac | 7 | #include "HMC5883L.h" |
| higedura | 0:813b917360ac | 8 | |
| higedura | 0:813b917360ac | 9 | Serial pc(USBTX, USBRX); |
| higedura | 0:813b917360ac | 10 | DigitalOut led1(LED1); |
| higedura | 0:813b917360ac | 11 | DigitalOut led2(LED2); |
| higedura | 0:813b917360ac | 12 | DigitalOut led3(LED3); |
| higedura | 0:813b917360ac | 13 | DigitalOut led4(LED4); |
| higedura | 0:813b917360ac | 14 | SDFileSystem sd(p5, p6, p7, p8, "sd"); // the pinout on the mbed Cool Components workshop board |
| higedura | 0:813b917360ac | 15 | ADXL345_I2C accelerometer(p9, p10); |
| higedura | 0:813b917360ac | 16 | ITG3200 gyro(p9, p10); |
| higedura | 0:813b917360ac | 17 | HMC5883L compass(p9, p10); |
| higedura | 0:813b917360ac | 18 | Serial xbee(p13, p14); |
| higedura | 0:813b917360ac | 19 | AnalogIn alt_in(p19); |
| higedura | 0:813b917360ac | 20 | DigitalIn stop(p20); |
| higedura | 0:813b917360ac | 21 | PwmOut ESC1(p21); |
| higedura | 0:813b917360ac | 22 | PwmOut ESC2(p22); |
| higedura | 0:813b917360ac | 23 | PwmOut ESC3(p23); |
| higedura | 0:813b917360ac | 24 | PwmOut ESC4(p24); |
| higedura | 0:813b917360ac | 25 | PwmOut ESC5(p25); |
| higedura | 0:813b917360ac | 26 | PwmOut ESC6(p26); |
| higedura | 0:813b917360ac | 27 | Serial navi(p28, p27); // tx, rx |
| higedura | 0:813b917360ac | 28 | DigitalOut navi_flag(p29); |
| higedura | 0:813b917360ac | 29 | |
| higedura | 0:813b917360ac | 30 | #define pi 3.14159265 |
| higedura | 0:813b917360ac | 31 | |
| higedura | 0:813b917360ac | 32 | #define N 3 // phi, the, psi |
| higedura | 0:813b917360ac | 33 | #define M 6 // The numbers of rotors |
| higedura | 0:813b917360ac | 34 | #define L 5 // The numbers of state quantities (optimal regulator) |
| higedura | 0:813b917360ac | 35 | #define chan 4 // The numbers of channels |
| higedura | 0:813b917360ac | 36 | #define chan_buf 17 |
| higedura | 0:813b917360ac | 37 | #define N_LWMA 100 |
| higedura | 0:813b917360ac | 38 | |
| higedura | 0:813b917360ac | 39 | int preparing_acc(); |
| higedura | 0:813b917360ac | 40 | double* calib(); |
| higedura | 0:813b917360ac | 41 | double* RK4( double, double[N], double[N] ); |
| higedura | 0:813b917360ac | 42 | double* func( double[N], double[N] ); |
| higedura | 0:813b917360ac | 43 | double* LWMA( double[N] ); |
| higedura | 0:813b917360ac | 44 | double* EKF_predict( double[N], double[N] ); |
| higedura | 0:813b917360ac | 45 | double* EKF_correct( double[N], double[N], double[N] ); |
| higedura | 0:813b917360ac | 46 | |
| higedura | 0:813b917360ac | 47 | // 0 1 2 |
| higedura | 0:813b917360ac | 48 | // [ accXn-1 accXn-2 ... ] 0 |
| higedura | 0:813b917360ac | 49 | // zBuf = [ accYn-1 accYn-2 ... ] 1 |
| higedura | 0:813b917360ac | 50 | // [ accZn-1 accZn-2 ... ] 2 |
| higedura | 0:813b917360ac | 51 | double z_buf[N][N_LWMA] = {0}; // For LWMA |
| higedura | 0:813b917360ac | 52 | |
| higedura | 0:813b917360ac | 53 | double P[N][N] = {{1,0,0},{0,1,0},{0,0,1}}; // For EKF |
| higedura | 0:813b917360ac | 54 | |
| higedura | 0:813b917360ac | 55 | int main(){ |
| higedura | 0:813b917360ac | 56 | |
| higedura | 0:813b917360ac | 57 | pc.baud(921600); |
| higedura | 0:813b917360ac | 58 | navi.baud(57600); |
| higedura | 0:813b917360ac | 59 | xbee.baud(57600); |
| higedura | 0:813b917360ac | 60 | |
| higedura | 0:813b917360ac | 61 | FILE *fp = fopen("/sd/sdtest.txt", "w"); |
| higedura | 0:813b917360ac | 62 | if(fp == NULL) { |
| higedura | 0:813b917360ac | 63 | error("Could not open file for write\n"); |
| higedura | 0:813b917360ac | 64 | while(1){ led2 = 1; wait(.5); led2 = 0; led3 = 0; wait(.5); led3 = 0; } |
| higedura | 0:813b917360ac | 65 | } |
| higedura | 0:813b917360ac | 66 | |
| higedura | 0:813b917360ac | 67 | int correct_period = 1000; |
| higedura | 0:813b917360ac | 68 | int navi_period = 10; |
| higedura | 0:813b917360ac | 69 | int xbee_period = 100; |
| higedura | 0:813b917360ac | 70 | int correct_loop = -2000; |
| higedura | 0:813b917360ac | 71 | int navi_loop = 10; |
| higedura | 0:813b917360ac | 72 | int xbee_loop = 100; |
| higedura | 0:813b917360ac | 73 | |
| higedura | 0:813b917360ac | 74 | double dt_wait = 0.0019; |
| higedura | 0:813b917360ac | 75 | //double dt_wait = 0.01; |
| higedura | 0:813b917360ac | 76 | double dt = 0.01; |
| higedura | 0:813b917360ac | 77 | double t = 0; |
| higedura | 0:813b917360ac | 78 | |
| higedura | 0:813b917360ac | 79 | int bit_acc[N] = {0}; // Buffer of the accelerometer |
| higedura | 0:813b917360ac | 80 | int get_mag[N] = {0}; // Buffer of the compass |
| higedura | 0:813b917360ac | 81 | |
| higedura | 0:813b917360ac | 82 | // Calibration routine |
| higedura | 0:813b917360ac | 83 | double calib_acc[N] = {0}; |
| higedura | 0:813b917360ac | 84 | double calib_gyro[N] = {0}; |
| higedura | 0:813b917360ac | 85 | double compass_basis_rad = 0; |
| higedura | 0:813b917360ac | 86 | |
| higedura | 0:813b917360ac | 87 | // Getting data |
| higedura | 0:813b917360ac | 88 | double acc[N] = {0, 0, 1}; |
| higedura | 0:813b917360ac | 89 | double acc0[N] = {0, 0, 1}; |
| higedura | 0:813b917360ac | 90 | double d_acc[N] = {0}; |
| higedura | 0:813b917360ac | 91 | double gyro_deg[N] = {0}; |
| higedura | 0:813b917360ac | 92 | double gyro_rad[N] = {0}; |
| higedura | 0:813b917360ac | 93 | int mag[N] = {0}; |
| higedura | 0:813b917360ac | 94 | |
| higedura | 0:813b917360ac | 95 | // Measurement algorithm |
| higedura | 0:813b917360ac | 96 | double angle_acc[2] = {0}; |
| higedura | 0:813b917360ac | 97 | double angle_deg[N] = {0}; |
| higedura | 0:813b917360ac | 98 | double angle_rad[N] = {0}; |
| higedura | 0:813b917360ac | 99 | double zLWMA[N] = {0}; |
| higedura | 0:813b917360ac | 100 | double compass_rad = 0; |
| higedura | 0:813b917360ac | 101 | double compass_deg = 0; |
| higedura | 0:813b917360ac | 102 | |
| higedura | 0:813b917360ac | 103 | // Gravity z |
| higedura | 0:813b917360ac | 104 | for( int i=0;i<N_LWMA;i++ ){ z_buf[2][i] = 1; } |
| higedura | 0:813b917360ac | 105 | |
| higedura | 0:813b917360ac | 106 | double* p_calib; |
| higedura | 0:813b917360ac | 107 | double* p_RK4; |
| higedura | 0:813b917360ac | 108 | double* p_EKF; |
| higedura | 0:813b917360ac | 109 | double* p_zLWMA; |
| higedura | 0:813b917360ac | 110 | |
| higedura | 0:813b917360ac | 111 | // Optimal regulator |
| higedura | 0:813b917360ac | 112 | double state[L]; |
| higedura | 0:813b917360ac | 113 | // Gain |
| higedura | 0:813b917360ac | 114 | double Kr[M][L] = {{ 0, 5.774, 0, 2.288, 0.408}, |
| higedura | 0:813b917360ac | 115 | { -5, 2.887,-1.982, 1.144,-0.408}, |
| higedura | 0:813b917360ac | 116 | { -5,-2.887,-1.982,-1.144, 0.408}, |
| higedura | 0:813b917360ac | 117 | { 0,-5.774, 0,-2.288,-0.408}, |
| higedura | 0:813b917360ac | 118 | { 5,-2.887, 1.982,-1.144, 0.408}, |
| higedura | 0:813b917360ac | 119 | { 5, 2.887, 1.982, 1.144,-0.408}}; |
| higedura | 0:813b917360ac | 120 | |
| higedura | 0:813b917360ac | 121 | double motor[M] = {0}; |
| higedura | 0:813b917360ac | 122 | double motor_attitude[M] = {0}; // Pulth width by the attitude control |
| higedura | 0:813b917360ac | 123 | double thrust[M]; |
| higedura | 0:813b917360ac | 124 | |
| higedura | 0:813b917360ac | 125 | // Hovering puls width |
| higedura | 0:813b917360ac | 126 | //double hover = 0.0015; |
| higedura | 0:813b917360ac | 127 | |
| higedura | 0:813b917360ac | 128 | // Coefficient of transfering Force to Puls |
| higedura | 0:813b917360ac | 129 | //double Cfp = 1/(1.4*pow(10.0, 6)); |
| higedura | 0:813b917360ac | 130 | double Cfp = 2*pow(10.0, -7); |
| higedura | 0:813b917360ac | 131 | // Coefficient of transfering altitude to Puls |
| higedura | 0:813b917360ac | 132 | double Calt = 1/(1.4*pow(10.0, 5)); |
| higedura | 0:813b917360ac | 133 | |
| higedura | 0:813b917360ac | 134 | // Keeping altitude |
| higedura | 0:813b917360ac | 135 | double alt; |
| higedura | 0:813b917360ac | 136 | double alt_r = 1000; // Altitude reference 1m |
| higedura | 0:813b917360ac | 137 | double alt_keep; |
| higedura | 0:813b917360ac | 138 | double Kp_alt = 1; |
| higedura | 0:813b917360ac | 139 | double sensor_h = 300; // Hight of sensor |
| higedura | 0:813b917360ac | 140 | |
| higedura | 0:813b917360ac | 141 | // Navigation |
| higedura | 0:813b917360ac | 142 | char command_buf1[chan_buf] = {0}; |
| higedura | 0:813b917360ac | 143 | int command_buf2[chan_buf] = {0}; |
| higedura | 0:813b917360ac | 144 | double command[chan] = {0}; |
| higedura | 0:813b917360ac | 145 | int set_navi_flag = 0; |
| higedura | 0:813b917360ac | 146 | |
| higedura | 0:813b917360ac | 147 | // *** Setting up sensors *** |
| higedura | 0:813b917360ac | 148 | //pc.printf("\r\n\r\nSetting up sensors\r\n"); |
| higedura | 0:813b917360ac | 149 | xbee.printf("\r\n\r\nSetting up sensors\r\n"); |
| higedura | 0:813b917360ac | 150 | |
| higedura | 0:813b917360ac | 151 | // These are here to test whether any of the initialization fails. It will print the failure. |
| higedura | 0:813b917360ac | 152 | if (accelerometer.setPowerControl(0x00)){ |
| higedura | 0:813b917360ac | 153 | //pc.printf("didn't intitialize power control\n\r"); |
| higedura | 0:813b917360ac | 154 | xbee.printf("didn't intitialize power control\n\r"); |
| higedura | 0:813b917360ac | 155 | return 0; |
| higedura | 0:813b917360ac | 156 | } |
| higedura | 0:813b917360ac | 157 | // Full resolution, +/-16g, 4mg/LSB. |
| higedura | 0:813b917360ac | 158 | wait(.1); |
| higedura | 0:813b917360ac | 159 | |
| higedura | 0:813b917360ac | 160 | if(accelerometer.setDataFormatControl(0x0B)){ |
| higedura | 0:813b917360ac | 161 | //pc.printf("didn't set data format\n\r"); |
| higedura | 0:813b917360ac | 162 | xbee.printf("didn't set data format\n\r"); |
| higedura | 0:813b917360ac | 163 | return 0; |
| higedura | 0:813b917360ac | 164 | } |
| higedura | 0:813b917360ac | 165 | wait(.1); |
| higedura | 0:813b917360ac | 166 | |
| higedura | 0:813b917360ac | 167 | // 3.2kHz data rate. |
| higedura | 0:813b917360ac | 168 | if(accelerometer.setDataRate(ADXL345_3200HZ)){ |
| higedura | 0:813b917360ac | 169 | //pc.printf("didn't set data rate\n\r"); |
| higedura | 0:813b917360ac | 170 | xbee.printf("didn't set data rate\n\r"); |
| higedura | 0:813b917360ac | 171 | return 0; |
| higedura | 0:813b917360ac | 172 | } |
| higedura | 0:813b917360ac | 173 | wait(.1); |
| higedura | 0:813b917360ac | 174 | |
| higedura | 0:813b917360ac | 175 | if(accelerometer.setPowerControl(MeasurementMode)) { |
| higedura | 0:813b917360ac | 176 | //pc.printf("didn't set the power control to measurement\n\r"); |
| higedura | 0:813b917360ac | 177 | xbee.printf("didn't set the power control to measurement\n\r"); |
| higedura | 0:813b917360ac | 178 | return 0; |
| higedura | 0:813b917360ac | 179 | } |
| higedura | 0:813b917360ac | 180 | wait(.1); |
| higedura | 0:813b917360ac | 181 | |
| higedura | 0:813b917360ac | 182 | gyro.setLpBandwidth(LPFBW_42HZ); |
| higedura | 0:813b917360ac | 183 | compass.setDefault(); |
| higedura | 0:813b917360ac | 184 | wait(.1); // Wait some time for all sensors (Need at least 5ms) |
| higedura | 0:813b917360ac | 185 | |
| higedura | 0:813b917360ac | 186 | // *** Setting up motors *** |
| higedura | 0:813b917360ac | 187 | //pc.printf("Setting up motors\r\n"); |
| higedura | 0:813b917360ac | 188 | xbee.printf("Setting up motors\r\n"); |
| higedura | 0:813b917360ac | 189 | ESC1.period(0.016044); ESC2.period(0.016044); ESC3.period(0.016044); ESC4.period(0.016044); ESC5.period(0.016044); ESC6.period(0.016044); |
| higedura | 0:813b917360ac | 190 | // pulsewidth 0.0011~0.00195 |
| higedura | 0:813b917360ac | 191 | ESC1.pulsewidth(0.001); ESC2.pulsewidth(0.001); ESC3.pulsewidth(0.001); ESC4.pulsewidth(0.001); ESC5.pulsewidth(0.001); ESC6.pulsewidth(0.001); |
| higedura | 0:813b917360ac | 192 | // Wait some time for ESC (about 10s) |
| higedura | 0:813b917360ac | 193 | wait(5); |
| higedura | 0:813b917360ac | 194 | |
| higedura | 0:813b917360ac | 195 | // *** Setting up navigator *** |
| higedura | 0:813b917360ac | 196 | //pc.printf("Setting up navigator\r\n"); |
| higedura | 0:813b917360ac | 197 | xbee.printf("Setting up navigator\r\n"); |
| higedura | 0:813b917360ac | 198 | while( set_navi_flag==0 ){ |
| higedura | 0:813b917360ac | 199 | navi_flag = 1; |
| higedura | 0:813b917360ac | 200 | for( int i=0;i<chan_buf;i++ ){ |
| higedura | 0:813b917360ac | 201 | navi.scanf("%c",&command_buf1[i]); |
| higedura | 0:813b917360ac | 202 | if(command_buf1[i]=='a'){ |
| higedura | 0:813b917360ac | 203 | set_navi_flag=1; |
| higedura | 0:813b917360ac | 204 | break; |
| higedura | 0:813b917360ac | 205 | } |
| higedura | 0:813b917360ac | 206 | } |
| higedura | 0:813b917360ac | 207 | navi_flag = 0; |
| higedura | 0:813b917360ac | 208 | wait(.1); |
| higedura | 0:813b917360ac | 209 | } |
| higedura | 0:813b917360ac | 210 | |
| higedura | 0:813b917360ac | 211 | // *** Calibration routine *** |
| higedura | 0:813b917360ac | 212 | p_calib = calib(); |
| higedura | 0:813b917360ac | 213 | for( int i=0;i<3;i++ ){ calib_acc[i] = *p_calib; p_calib = p_calib+1; } |
| higedura | 0:813b917360ac | 214 | for( int i=3;i<6;i++ ){ calib_gyro[i-3] = *p_calib; p_calib = p_calib+1; } |
| higedura | 0:813b917360ac | 215 | compass_basis_rad = *p_calib; |
| higedura | 0:813b917360ac | 216 | |
| higedura | 0:813b917360ac | 217 | led1 = 1; led2 = 1; led3 = 1; led4 = 1; |
| higedura | 0:813b917360ac | 218 | |
| higedura | 0:813b917360ac | 219 | //pc.printf("Starting IMU unit\n\r"); |
| higedura | 0:813b917360ac | 220 | xbee.printf("Starting IMU unit\n\r"); |
| higedura | 0:813b917360ac | 221 | //pc.printf(" time phi the P Q R accX accY accZ alt T1 T2 T3 T4 T5 T6 thro roll pitch yaw\n\r"); |
| higedura | 0:813b917360ac | 222 | xbee.printf(" time phi the P Q R accX accY accZ alt T1 T2 T3 T4 T5 T6 thro roll pitch yaw\n\r"); |
| higedura | 0:813b917360ac | 223 | |
| higedura | 0:813b917360ac | 224 | //while(1){ |
| higedura | 0:813b917360ac | 225 | while( stop==0 ){ |
| higedura | 0:813b917360ac | 226 | |
| higedura | 0:813b917360ac | 227 | // Navigation |
| higedura | 0:813b917360ac | 228 | if( navi_loop>=navi_period ){ |
| higedura | 0:813b917360ac | 229 | |
| higedura | 0:813b917360ac | 230 | navi_flag = 1; |
| higedura | 0:813b917360ac | 231 | |
| higedura | 0:813b917360ac | 232 | for( int i=0;i<chan_buf;i++ ){ navi.scanf("%c",&command_buf1[i]); } |
| higedura | 0:813b917360ac | 233 | for( int i=0;i<chan_buf;i++ ){ |
| higedura | 0:813b917360ac | 234 | command_buf2[i] = (int)command_buf1[i]-48; |
| higedura | 0:813b917360ac | 235 | if( command_buf2[i]==-16 ){ command_buf2[i]=0; } |
| higedura | 0:813b917360ac | 236 | } |
| higedura | 0:813b917360ac | 237 | for( int i=0;i<chan;i++ ){ command[i] = 0.0001*command_buf2[1+i*4]+0.00001*command_buf2[2+i*4]+0.000001*command_buf2[3+i*4]+0.001; } |
| higedura | 0:813b917360ac | 238 | |
| higedura | 0:813b917360ac | 239 | navi_flag = 0; |
| higedura | 0:813b917360ac | 240 | navi_loop = 0; |
| higedura | 0:813b917360ac | 241 | |
| higedura | 0:813b917360ac | 242 | } |
| higedura | 0:813b917360ac | 243 | navi_loop++; |
| higedura | 0:813b917360ac | 244 | //command[0] = 0.0016; // Shiyougo comentout surukoto!! |
| higedura | 0:813b917360ac | 245 | |
| higedura | 0:813b917360ac | 246 | // Updating accelerometer and compass |
| higedura | 0:813b917360ac | 247 | accelerometer.getOutput(bit_acc); |
| higedura | 0:813b917360ac | 248 | compass.readData(get_mag); |
| higedura | 0:813b917360ac | 249 | alt = alt_in/0.000074; |
| higedura | 0:813b917360ac | 250 | |
| higedura | 0:813b917360ac | 251 | for( int i=0;i<N;i++ ){ acc0[i] = acc[i]; } |
| higedura | 0:813b917360ac | 252 | |
| higedura | 0:813b917360ac | 253 | // Transfering units and Coordinate transform |
| higedura | 0:813b917360ac | 254 | acc[0] = (((int16_t)bit_acc[0]-calib_acc[0])-((int16_t)bit_acc[1]-calib_acc[1]))*0.004/sqrt(2.0); |
| higedura | 0:813b917360ac | 255 | acc[1] = (((int16_t)bit_acc[0]-calib_acc[0])+((int16_t)bit_acc[1]-calib_acc[1]))*0.004/sqrt(2.0); |
| higedura | 0:813b917360ac | 256 | acc[2] = ((int16_t)bit_acc[2]-calib_acc[2])*0.004+1; |
| higedura | 0:813b917360ac | 257 | |
| higedura | 0:813b917360ac | 258 | gyro_deg[0] = ((gyro.getGyroX()-calib_gyro[0])-(gyro.getGyroY()-calib_gyro[1]))/14.375/sqrt(2.0); |
| higedura | 0:813b917360ac | 259 | gyro_deg[1] = (-(gyro.getGyroX()-calib_gyro[0])-(gyro.getGyroY()-calib_gyro[1]))/14.375/sqrt(2.0); // Modify the differencial of the sensor axis |
| higedura | 0:813b917360ac | 260 | gyro_deg[2] = -(gyro.getGyroZ()-calib_gyro[2])/14.375; |
| higedura | 0:813b917360ac | 261 | |
| higedura | 0:813b917360ac | 262 | for( int i=0;i<N;i++ ){ mag[0] = (int16_t)get_mag[0]; } |
| higedura | 0:813b917360ac | 263 | |
| higedura | 0:813b917360ac | 264 | // Low pass filter for acc |
| higedura | 0:813b917360ac | 265 | //if( -0.05<acc[0] && acc[0]<0.05 ){ acc[0]=0; } |
| higedura | 0:813b917360ac | 266 | //if( -0.05<acc[1] && acc[1]<0.05 ){ acc[1]=0; } |
| higedura | 0:813b917360ac | 267 | //if( 0.95<acc[2] && acc[2]<1.05 ){ acc[2]=0; } |
| higedura | 0:813b917360ac | 268 | // Limitter for acc |
| higedura | 0:813b917360ac | 269 | if( acc[0]<-1 ){ acc[0]=-1; } if( 1<acc[0] ){ acc[0]=1; } |
| higedura | 0:813b917360ac | 270 | if( acc[1]<-1 ){ acc[1]=-1; } if( 1<acc[1] ){ acc[1]=1; } |
| higedura | 0:813b917360ac | 271 | if( acc[2]<-0 ){ acc[2]=-0; } if( 2<acc[2] ){ acc[2]=2; } |
| higedura | 0:813b917360ac | 272 | |
| higedura | 0:813b917360ac | 273 | for( int i=0;i<N;i++ ){ d_acc[i] = acc[i]-acc0[i]; } |
| higedura | 0:813b917360ac | 274 | if( abs(d_acc[0])>=0.5 && abs(acc[0])>=0.3 ){ acc[0] = acc0[0]; } |
| higedura | 0:813b917360ac | 275 | if( abs(d_acc[1])>=0.5 && abs(acc[1])>=0.3 ){ acc[1] = acc0[1]; } |
| higedura | 0:813b917360ac | 276 | if( abs(d_acc[2])>=0.5 ){ acc[2] = acc0[2]; } |
| higedura | 0:813b917360ac | 277 | |
| higedura | 0:813b917360ac | 278 | // Low pass filter for gyro |
| higedura | 0:813b917360ac | 279 | //for( int i=0;i<N;i++ ){if( -0.5<gyro_deg[i] && gyro_deg[i]<0.5 ){ gyro_deg[i] = 0; }} |
| higedura | 0:813b917360ac | 280 | // Limitter for gyro |
| higedura | 0:813b917360ac | 281 | //for( int i=0;i<N;i++ ){if( gyro_deg[i]<-90 ){ gyro_deg[i]=-90; }} for( int i=0;i<N;i++ ){if( 90<gyro_deg[i] ){ gyro_deg[i]=90; }} |
| higedura | 0:813b917360ac | 282 | |
| higedura | 0:813b917360ac | 283 | for( int i=0;i<N;i++ ){ gyro_rad[i] = gyro_deg[i]*pi/180; } |
| higedura | 0:813b917360ac | 284 | /* |
| higedura | 0:813b917360ac | 285 | // Compass yaw |
| higedura | 0:813b917360ac | 286 | compass_rad = (double)mag[1]/mag[0]; |
| higedura | 0:813b917360ac | 287 | compass_rad = atan(compass_rad); |
| higedura | 0:813b917360ac | 288 | //compass_rad = compass_rad-compass_basis_rad; |
| higedura | 0:813b917360ac | 289 | compass_deg = compass_rad*180/pi; |
| higedura | 0:813b917360ac | 290 | */ |
| higedura | 0:813b917360ac | 291 | // LWMA (Observation) |
| higedura | 0:813b917360ac | 292 | p_zLWMA = LWMA(acc); |
| higedura | 0:813b917360ac | 293 | for( int i=0;i<N;i++ ){ zLWMA[i] = *p_zLWMA; p_zLWMA = p_zLWMA+1; } |
| higedura | 0:813b917360ac | 294 | // LWMA angle |
| higedura | 0:813b917360ac | 295 | //angle_acc[0] = asin(zLWMA[1])*180/pi; |
| higedura | 0:813b917360ac | 296 | //angle_acc[1] = asin(zLWMA[0])*180/pi; |
| higedura | 0:813b917360ac | 297 | |
| higedura | 0:813b917360ac | 298 | // RK4 (Prediction) |
| higedura | 0:813b917360ac | 299 | p_RK4 = RK4(dt,angle_rad,gyro_rad); |
| higedura | 0:813b917360ac | 300 | for( int i=0;i<N;i++ ){ angle_rad[i] = *p_RK4; p_RK4 = p_RK4+1; } |
| higedura | 0:813b917360ac | 301 | |
| higedura | 0:813b917360ac | 302 | // EKF (Correction) |
| higedura | 0:813b917360ac | 303 | EKF_predict(angle_rad,gyro_rad); |
| higedura | 0:813b917360ac | 304 | if ( correct_loop>=correct_period ){ |
| higedura | 0:813b917360ac | 305 | p_EKF = EKF_correct(angle_rad,gyro_rad,zLWMA); |
| higedura | 0:813b917360ac | 306 | for ( int i=0;i<N;i++ ){ angle_rad[i] = *p_EKF; p_EKF = p_EKF+1; } |
| higedura | 0:813b917360ac | 307 | correct_loop = 0; |
| higedura | 0:813b917360ac | 308 | } |
| higedura | 0:813b917360ac | 309 | correct_loop++; |
| higedura | 0:813b917360ac | 310 | |
| higedura | 0:813b917360ac | 311 | // Caribrating EKF |
| higedura | 0:813b917360ac | 312 | for( int i=0;i<N;i++ ){ angle_deg[i] = angle_rad[i]*180/pi; } |
| higedura | 0:813b917360ac | 313 | |
| higedura | 0:813b917360ac | 314 | alt_keep = -Kp_alt*(alt-alt_r-sensor_h)*Calt; |
| higedura | 0:813b917360ac | 315 | //alt_keep = 0; |
| higedura | 0:813b917360ac | 316 | |
| higedura | 0:813b917360ac | 317 | state[0] = angle_deg[0]; state[1] = angle_deg[1]; |
| higedura | 0:813b917360ac | 318 | state[2] = gyro_deg[0]; state[3] = gyro_deg[1]; state[4] = gyro_deg[2]; |
| higedura | 0:813b917360ac | 319 | |
| higedura | 0:813b917360ac | 320 | for( int i=0;i<M;i++ ){ |
| higedura | 0:813b917360ac | 321 | motor_attitude[i] = 0; |
| higedura | 0:813b917360ac | 322 | for( int j=0;j<L;j++ ){ motor_attitude[i] += -Kr[i][j]*state[j]; } |
| higedura | 0:813b917360ac | 323 | } |
| higedura | 0:813b917360ac | 324 | |
| higedura | 0:813b917360ac | 325 | for( int i=0;i<M;i++ ){ motor[i] = command[0]+motor_attitude[i]*Cfp; } |
| higedura | 0:813b917360ac | 326 | |
| higedura | 0:813b917360ac | 327 | // pulsewidth 0.0011~0.0019 (0.001~0.002?) |
| higedura | 0:813b917360ac | 328 | for( int i=0;i<M;i++ ){ |
| higedura | 0:813b917360ac | 329 | if( motor[i]>0.00195 ){ motor[i]=0.00195; } |
| higedura | 0:813b917360ac | 330 | if( motor[i]<0.0011 ){ motor[i]=0.00105; } |
| higedura | 0:813b917360ac | 331 | if( command[0]<0.0011 ){ motor[i]=0.00105; } |
| higedura | 0:813b917360ac | 332 | } |
| higedura | 0:813b917360ac | 333 | |
| higedura | 0:813b917360ac | 334 | ESC1.pulsewidth(motor[0]); ESC2.pulsewidth(motor[1]); ESC3.pulsewidth(motor[2]); ESC4.pulsewidth(motor[3]); ESC5.pulsewidth(motor[4]); ESC6.pulsewidth(motor[5]); |
| higedura | 0:813b917360ac | 335 | |
| higedura | 0:813b917360ac | 336 | for( int i=0;i<M;i++ ){ thrust[i] = 630000.0*motor[i]*motor[i]-700.0*motor[i]; } |
| higedura | 0:813b917360ac | 337 | |
| higedura | 0:813b917360ac | 338 | //pc.printf("%7.2f, %7.3f, %7.3f, %7.3f, %7.1f, %7.1f, %7.1f, %5d, %5d, %5d\n\r", t, acc[0], acc[1], acc[2], gyro_deg[0], gyro_deg[1], gyro_deg[2], mag[0], mag[1], mag[2]); |
| higedura | 0:813b917360ac | 339 | //pc.printf("%7.2f, %7.1f, %7.1f, %7.1f, %7.1f, %7.1f, %7.3f, %7.3f, %7.3f, %7.1f, %6.3f, %6.3f, %6.3f, %6.3f, %6.3f, %6.3f, %6.3f, %6.3f, %6.3f, %6.3f\n\r", t, state[0], state[1], state[2], state[3], state[4], zLWMA[0], zLWMA[1], zLWMA[2], alt, motor[0]*1000, motor[1]*1000, motor[2]*1000, motor[3]*1000, motor[4]*1000, motor[5]*1000, command[0]*1000, command[1]*1000, command[2]*1000, command[3]*1000); |
| higedura | 0:813b917360ac | 340 | fprintf(fp,"%7.2f, %7.1f, %7.1f, %7.1f, %7.1f, %7.1f, %7.3f, %7.3f, %7.3f, %7.3f, %7.3f, %7.3f, %7.1f, %6.3f, %6.3f, %6.3f, %6.3f, %6.3f, %6.3f, %6.3f, %6.3f, %6.3f, %6.3f\n", t, state[0], state[1], state[2], state[3], state[4], acc[0], acc[1], acc[2], zLWMA[0], zLWMA[1], zLWMA[2], alt, motor[0]*1000, motor[1]*1000, motor[2]*1000, motor[3]*1000, motor[4]*1000, motor[5]*1000, command[0]*1000, command[1]*1000, command[2]*1000, command[3]*1000); |
| higedura | 0:813b917360ac | 341 | if( xbee_loop>=xbee_period ){ |
| higedura | 0:813b917360ac | 342 | xbee.printf("%7.2f, %7.1f, %7.1f, %7.1f, %7.1f, %7.1f, %7.3f, %7.3f, %7.3f, %7.3f, %7.3f, %7.3f, %7.1f, %6.3f, %6.3f, %6.3f, %6.3f, %6.3f, %6.3f, %6.3f, %6.3f, %6.3f, %6.3f\n\r", t, state[0], state[1], state[2], state[3], state[4], acc[0], acc[1], acc[2], zLWMA[0], zLWMA[1], zLWMA[2], alt, motor[0]*1000, motor[1]*1000, motor[2]*1000, motor[3]*1000, motor[4]*1000, motor[5]*1000, command[0]*1000, command[1]*1000, command[2]*1000, command[3]*1000); |
| higedura | 0:813b917360ac | 343 | xbee_loop = 0; |
| higedura | 0:813b917360ac | 344 | } |
| higedura | 0:813b917360ac | 345 | xbee_loop++; |
| higedura | 0:813b917360ac | 346 | wait(dt_wait); |
| higedura | 0:813b917360ac | 347 | t += dt; |
| higedura | 0:813b917360ac | 348 | |
| higedura | 0:813b917360ac | 349 | } |
| higedura | 0:813b917360ac | 350 | |
| higedura | 0:813b917360ac | 351 | // pulsewidth 0.0011~0.00195 |
| higedura | 0:813b917360ac | 352 | ESC1.pulsewidth(0.001); ESC2.pulsewidth(0.001); ESC3.pulsewidth(0.001); ESC4.pulsewidth(0.001); ESC5.pulsewidth(0.001); ESC6.pulsewidth(0.001); |
| higedura | 0:813b917360ac | 353 | led1 = 0; led2 = 0; led3 = 0; led4 = 0; |
| higedura | 0:813b917360ac | 354 | fclose(fp); |
| higedura | 0:813b917360ac | 355 | |
| higedura | 0:813b917360ac | 356 | } |
| higedura | 0:813b917360ac | 357 | |
| higedura | 0:813b917360ac | 358 | double* EKF_predict( double y[N], double x[N] ){ |
| higedura | 0:813b917360ac | 359 | // x = F * x; |
| higedura | 0:813b917360ac | 360 | // P = F * P * F' + G * Q * G'; |
| higedura | 0:813b917360ac | 361 | |
| higedura | 0:813b917360ac | 362 | //double Q[N][N] = { {0.1, 0, 0}, {0, 0.1, 0}, {0, 0, 0.1} }; |
| higedura | 0:813b917360ac | 363 | double Q[N][N] = { {0.00263, 0, 0}, {0, 0.00373, 0}, {0, 0, 0.00509} }; |
| higedura | 0:813b917360ac | 364 | |
| higedura | 0:813b917360ac | 365 | double Fjac[N][N] = {{x[1]*cos(y[0])*tan(y[1])-x[2]*sin(y[0])*tan(y[1]), x[1]*sin(y[0])/(cos(y[1])*cos(y[1]))+x[2]*cos(y[0])/(cos(y[1])*cos(y[1])), 0}, {-x[1]*sin(y[0])-x[2]*cos(y[0]), 0, 0}, {x[1]*cos(y[0])/cos(y[1])-x[2]*sin(y[0])/cos(y[1]), x[1]*sin(y[0])*sin(y[1])/(cos(y[1])*cos(y[1]))+x[2]*cos(y[0])*sin(y[1])/(cos(y[1])*cos(y[1])), 0}}; |
| higedura | 0:813b917360ac | 366 | double Fjac_t[N][N] = {{x[1]*cos(y[0])*tan(y[1])-x[2]*sin(y[0])*tan(y[1]), -x[1]*sin(y[0])-x[2]*cos(y[0]), x[1]*cos(y[0])/cos(y[1])-x[2]*sin(y[0])/cos(y[1])}, {x[1]*sin(y[0])/(cos(y[1])*cos(y[1]))+x[2]*cos(y[0])/(cos(y[1])*cos(y[1])), 0, x[1]*sin(y[0])*sin(y[1])/(cos(y[1])*cos(y[1]))+x[2]*cos(y[0])*sin(y[1])/(cos(y[1])*cos(y[1]))}, {0, 0, 0}}; |
| higedura | 0:813b917360ac | 367 | double Gjac[N][N] = {{1, sin(y[0])*tan(y[1]), cos(y[0])*tan(y[1])}, {0, cos(y[0]), -sin(y[0])}, {0, sin(y[0])/cos(y[1]), cos(y[0])/cos(y[1])}}; |
| higedura | 0:813b917360ac | 368 | double Gjac_t[N][N] = {{1, 0, 0}, {sin(y[0])*tan(y[1]), cos(y[0]), sin(y[0])/cos(y[1])}, {cos(y[0])*tan(y[1]), -sin(y[0]), cos(y[0])/cos(y[1])}}; |
| higedura | 0:813b917360ac | 369 | |
| higedura | 0:813b917360ac | 370 | double FPF[N][N] = {0}; |
| higedura | 0:813b917360ac | 371 | double GQG[N][N] = {0}; |
| higedura | 0:813b917360ac | 372 | |
| higedura | 0:813b917360ac | 373 | double FP[N][N] = {0}; |
| higedura | 0:813b917360ac | 374 | double GQ[N][N] = {0}; |
| higedura | 0:813b917360ac | 375 | |
| higedura | 0:813b917360ac | 376 | for( int i=0;i<N;i++ ){ |
| higedura | 0:813b917360ac | 377 | for( int j=0;j<N;j++ ){ |
| higedura | 0:813b917360ac | 378 | for( int k=0;k<N;k++ ){ |
| higedura | 0:813b917360ac | 379 | FP[i][j] += Fjac[i][k]*P[k][j]; |
| higedura | 0:813b917360ac | 380 | GQ[i][j] += Gjac[i][k]*Q[k][j]; |
| higedura | 0:813b917360ac | 381 | |
| higedura | 0:813b917360ac | 382 | } |
| higedura | 0:813b917360ac | 383 | } |
| higedura | 0:813b917360ac | 384 | } |
| higedura | 0:813b917360ac | 385 | |
| higedura | 0:813b917360ac | 386 | for( int i=0;i<N;i++ ){ |
| higedura | 0:813b917360ac | 387 | for( int j=0;j<N;j++ ){ |
| higedura | 0:813b917360ac | 388 | for( int k=0;k<N;k++ ){ |
| higedura | 0:813b917360ac | 389 | FPF[i][j] += FP[i][k]*Fjac_t[k][j]; |
| higedura | 0:813b917360ac | 390 | GQG[i][j] += GQ[i][k]*Gjac_t[k][j]; |
| higedura | 0:813b917360ac | 391 | } |
| higedura | 0:813b917360ac | 392 | } |
| higedura | 0:813b917360ac | 393 | } |
| higedura | 0:813b917360ac | 394 | for( int i=0;i<N;i++ ){ |
| higedura | 0:813b917360ac | 395 | for( int j=0;j<N;j++ ){ |
| higedura | 0:813b917360ac | 396 | P[i][j] = FPF[i][j]+GQG[i][j]; |
| higedura | 0:813b917360ac | 397 | } |
| higedura | 0:813b917360ac | 398 | } |
| higedura | 0:813b917360ac | 399 | |
| higedura | 0:813b917360ac | 400 | return 0; |
| higedura | 0:813b917360ac | 401 | |
| higedura | 0:813b917360ac | 402 | } |
| higedura | 0:813b917360ac | 403 | |
| higedura | 0:813b917360ac | 404 | double* EKF_correct( double y[N], double x[N], double z[N] ){ |
| higedura | 0:813b917360ac | 405 | // K = P * H' / ( H * P * H' + R ) |
| higedura | 0:813b917360ac | 406 | // x = x + K * ( yobs(t) - H * x ) |
| higedura | 0:813b917360ac | 407 | // P = P - K * H * P |
| higedura | 0:813b917360ac | 408 | |
| higedura | 0:813b917360ac | 409 | //double R[N][N] = { {0.15, 0, 0}, {0, 0.15, 0}, {0, 0, 0.15} }; |
| higedura | 0:813b917360ac | 410 | double R[N][N] = { {0.00015, 0, 0}, {0, 0.00016, 0}, {0, 0, 0.00032} }; |
| higedura | 0:813b917360ac | 411 | |
| higedura | 0:813b917360ac | 412 | double Hjac[N][N] = {{0, cos(y[1]), 0}, {cos(y[0]), 0, 0}, {-sin(y[0])*cos(y[1]), -cos(y[0])*sin(y[1]), 0}}; |
| higedura | 0:813b917360ac | 413 | double Hjac_t[N][N] = {{0, cos(y[0]), -sin(y[0])*cos(y[1])}, {cos(y[1]), 0, -cos(y[0])*sin(y[1])}, {0, 0, 0}}; |
| higedura | 0:813b917360ac | 414 | double K[N][N] = {0}; // Kalman gain |
| higedura | 0:813b917360ac | 415 | double K_deno[N][N] = {0}; // Denominator of the kalman gain |
| higedura | 0:813b917360ac | 416 | double det_K_deno_inv = 0; |
| higedura | 0:813b917360ac | 417 | double K_deno_inv[N][N] = {0}; |
| higedura | 0:813b917360ac | 418 | double HPH[N][N] = {0}; |
| higedura | 0:813b917360ac | 419 | double HP[N][N] = {0}; |
| higedura | 0:813b917360ac | 420 | double PH[N][N] = {0}; |
| higedura | 0:813b917360ac | 421 | double KHP[N][N] = {0}; |
| higedura | 0:813b917360ac | 422 | |
| higedura | 0:813b917360ac | 423 | double Hx[N] = {0}; |
| higedura | 0:813b917360ac | 424 | double z_Hx[N] = {0}; |
| higedura | 0:813b917360ac | 425 | double Kz_Hx[N] = {0}; |
| higedura | 0:813b917360ac | 426 | |
| higedura | 0:813b917360ac | 427 | double* py = y; |
| higedura | 0:813b917360ac | 428 | |
| higedura | 0:813b917360ac | 429 | // Kalman gain |
| higedura | 0:813b917360ac | 430 | for( int i=0;i<N;i++ ){ |
| higedura | 0:813b917360ac | 431 | for( int j=0;j<N;j++ ){ |
| higedura | 0:813b917360ac | 432 | for( int k=0;k<N;k++ ){ |
| higedura | 0:813b917360ac | 433 | HP[i][j] += Hjac[i][k]*P[k][j]; |
| higedura | 0:813b917360ac | 434 | PH[i][j] += P[i][k]*Hjac_t[k][j]; |
| higedura | 0:813b917360ac | 435 | } |
| higedura | 0:813b917360ac | 436 | } |
| higedura | 0:813b917360ac | 437 | } |
| higedura | 0:813b917360ac | 438 | for( int i=0;i<N;i++ ){ |
| higedura | 0:813b917360ac | 439 | for( int j=0;j<N;j++ ){ |
| higedura | 0:813b917360ac | 440 | for( int k=0;k<N;k++ ){ |
| higedura | 0:813b917360ac | 441 | HPH[i][j] += HP[i][k]*Hjac_t[k][j]; |
| higedura | 0:813b917360ac | 442 | } |
| higedura | 0:813b917360ac | 443 | } |
| higedura | 0:813b917360ac | 444 | } |
| higedura | 0:813b917360ac | 445 | for( int i=0;i<N;i++ ){ |
| higedura | 0:813b917360ac | 446 | for( int j=0;j<N;j++ ){ |
| higedura | 0:813b917360ac | 447 | K_deno[i][j] = HPH[i][j]+R[i][j]; |
| higedura | 0:813b917360ac | 448 | } |
| higedura | 0:813b917360ac | 449 | } |
| higedura | 0:813b917360ac | 450 | // Inverce matrix |
| higedura | 0:813b917360ac | 451 | det_K_deno_inv = K_deno[0][0]*K_deno[1][1]*K_deno[2][2]+K_deno[1][0]*K_deno[2][1]*K_deno[0][2]+K_deno[2][0]*K_deno[0][1]*K_deno[1][2]-K_deno[0][0]*K_deno[2][1]*K_deno[1][2]-K_deno[2][0]*K_deno[1][1]*K_deno[0][2]-K_deno[1][0]*K_deno[0][1]*K_deno[2][2]; |
| higedura | 0:813b917360ac | 452 | K_deno_inv[0][0] = (K_deno[1][1]*K_deno[2][2]-K_deno[1][2]*K_deno[2][1])/det_K_deno_inv; |
| higedura | 0:813b917360ac | 453 | K_deno_inv[0][1] = (K_deno[0][2]*K_deno[2][1]-K_deno[0][1]*K_deno[2][2])/det_K_deno_inv; |
| higedura | 0:813b917360ac | 454 | K_deno_inv[0][2] = (K_deno[0][1]*K_deno[1][2]-K_deno[0][2]*K_deno[1][1])/det_K_deno_inv; |
| higedura | 0:813b917360ac | 455 | K_deno_inv[1][0] = (K_deno[1][2]*K_deno[2][0]-K_deno[1][0]*K_deno[2][2])/det_K_deno_inv; |
| higedura | 0:813b917360ac | 456 | K_deno_inv[1][1] = (K_deno[0][0]*K_deno[2][2]-K_deno[0][2]*K_deno[2][0])/det_K_deno_inv; |
| higedura | 0:813b917360ac | 457 | K_deno_inv[1][2] = (K_deno[0][2]*K_deno[1][0]-K_deno[0][0]*K_deno[1][2])/det_K_deno_inv; |
| higedura | 0:813b917360ac | 458 | K_deno_inv[2][0] = (K_deno[1][0]*K_deno[2][1]-K_deno[1][1]*K_deno[2][0])/det_K_deno_inv; |
| higedura | 0:813b917360ac | 459 | K_deno_inv[2][1] = (K_deno[0][1]*K_deno[2][0]-K_deno[0][0]*K_deno[2][1])/det_K_deno_inv; |
| higedura | 0:813b917360ac | 460 | K_deno_inv[2][2] = (K_deno[0][0]*K_deno[1][1]-K_deno[0][1]*K_deno[1][0])/det_K_deno_inv; |
| higedura | 0:813b917360ac | 461 | |
| higedura | 0:813b917360ac | 462 | for( int i=0;i<N;i++ ){ |
| higedura | 0:813b917360ac | 463 | for( int j=0;j<N;j++ ){ |
| higedura | 0:813b917360ac | 464 | for( int k=0;k<N;k++ ){ |
| higedura | 0:813b917360ac | 465 | K[i][j] += PH[i][k]*K_deno_inv[k][j]; |
| higedura | 0:813b917360ac | 466 | } |
| higedura | 0:813b917360ac | 467 | } |
| higedura | 0:813b917360ac | 468 | } |
| higedura | 0:813b917360ac | 469 | |
| higedura | 0:813b917360ac | 470 | // Filtering |
| higedura | 0:813b917360ac | 471 | for( int i=0;i<N;i++ ){ |
| higedura | 0:813b917360ac | 472 | for( int j=0;j<N;j++ ){ |
| higedura | 0:813b917360ac | 473 | Hx[i] += Hjac[i][j]*y[j]; |
| higedura | 0:813b917360ac | 474 | } |
| higedura | 0:813b917360ac | 475 | } |
| higedura | 0:813b917360ac | 476 | for( int i=0;i<N;i++ ){ |
| higedura | 0:813b917360ac | 477 | z_Hx[i] = z[i]-Hx[i]; |
| higedura | 0:813b917360ac | 478 | } |
| higedura | 0:813b917360ac | 479 | for( int i=0;i<N;i++ ){ |
| higedura | 0:813b917360ac | 480 | for( int j=0;j<N;j++ ){ |
| higedura | 0:813b917360ac | 481 | Kz_Hx[i] += K[i][j]*z_Hx[j]; |
| higedura | 0:813b917360ac | 482 | } |
| higedura | 0:813b917360ac | 483 | } |
| higedura | 0:813b917360ac | 484 | for( int i=0;i<N;i++ ){ |
| higedura | 0:813b917360ac | 485 | y[i] = y[i]+Kz_Hx[i]; |
| higedura | 0:813b917360ac | 486 | } |
| higedura | 0:813b917360ac | 487 | |
| higedura | 0:813b917360ac | 488 | // Covarience |
| higedura | 0:813b917360ac | 489 | for( int i=0;i<N;i++ ){ |
| higedura | 0:813b917360ac | 490 | for( int j=0;j<N;j++ ){ |
| higedura | 0:813b917360ac | 491 | for( int k=0;k<N;k++ ){ |
| higedura | 0:813b917360ac | 492 | KHP[i][j] += K[i][k]*HP[k][j]; |
| higedura | 0:813b917360ac | 493 | } |
| higedura | 0:813b917360ac | 494 | } |
| higedura | 0:813b917360ac | 495 | } |
| higedura | 0:813b917360ac | 496 | for( int i=0;i<N;i++ ){ |
| higedura | 0:813b917360ac | 497 | for( int j=0;j<N;j++ ){ |
| higedura | 0:813b917360ac | 498 | P[i][j] = P[i][j]-KHP[i][j]; |
| higedura | 0:813b917360ac | 499 | } |
| higedura | 0:813b917360ac | 500 | } |
| higedura | 0:813b917360ac | 501 | |
| higedura | 0:813b917360ac | 502 | return py; |
| higedura | 0:813b917360ac | 503 | |
| higedura | 0:813b917360ac | 504 | } |
| higedura | 0:813b917360ac | 505 | |
| higedura | 0:813b917360ac | 506 | double* LWMA( double z[N] ){ |
| higedura | 0:813b917360ac | 507 | |
| higedura | 0:813b917360ac | 508 | double zLWMA[N] = {0}; |
| higedura | 0:813b917360ac | 509 | double zLWMA_num[N] = {0}; |
| higedura | 0:813b917360ac | 510 | double zLWMA_den = 0; |
| higedura | 0:813b917360ac | 511 | |
| higedura | 0:813b917360ac | 512 | double* pzLWMA = zLWMA; |
| higedura | 0:813b917360ac | 513 | |
| higedura | 0:813b917360ac | 514 | for( int i=1;i<N_LWMA;i++ ){ |
| higedura | 0:813b917360ac | 515 | for( int j=0;j<N;j++ ){ |
| higedura | 0:813b917360ac | 516 | z_buf[j][N_LWMA-i] = z_buf[j][N_LWMA-i-1]; |
| higedura | 0:813b917360ac | 517 | } |
| higedura | 0:813b917360ac | 518 | } |
| higedura | 0:813b917360ac | 519 | for( int i=0;i<N;i++ ){ |
| higedura | 0:813b917360ac | 520 | z_buf[i][0] = z[i]; |
| higedura | 0:813b917360ac | 521 | } |
| higedura | 0:813b917360ac | 522 | |
| higedura | 0:813b917360ac | 523 | for( int i=0;i<N_LWMA;i++ ){ |
| higedura | 0:813b917360ac | 524 | for( int j=0;j<N;j++ ){ |
| higedura | 0:813b917360ac | 525 | zLWMA_num[j] += (N_LWMA-i)*z_buf[j][i]; |
| higedura | 0:813b917360ac | 526 | } |
| higedura | 0:813b917360ac | 527 | zLWMA_den += N_LWMA-i; |
| higedura | 0:813b917360ac | 528 | } |
| higedura | 0:813b917360ac | 529 | for( int i=0;i<N;i++ ){ |
| higedura | 0:813b917360ac | 530 | zLWMA[i] = zLWMA_num[i]/zLWMA_den; |
| higedura | 0:813b917360ac | 531 | } |
| higedura | 0:813b917360ac | 532 | |
| higedura | 0:813b917360ac | 533 | return pzLWMA; |
| higedura | 0:813b917360ac | 534 | |
| higedura | 0:813b917360ac | 535 | } |
| higedura | 0:813b917360ac | 536 | |
| higedura | 0:813b917360ac | 537 | double* RK4( double dt, double y[N], double x[N] ){ |
| higedura | 0:813b917360ac | 538 | |
| higedura | 0:813b917360ac | 539 | double yBuf[N] = {0}; |
| higedura | 0:813b917360ac | 540 | double k[N][4] = {0}; |
| higedura | 0:813b917360ac | 541 | |
| higedura | 0:813b917360ac | 542 | double* p_y = y; |
| higedura | 0:813b917360ac | 543 | |
| higedura | 0:813b917360ac | 544 | double* pk1; |
| higedura | 0:813b917360ac | 545 | double* pk2; |
| higedura | 0:813b917360ac | 546 | double* pk3; |
| higedura | 0:813b917360ac | 547 | double* pk4; |
| higedura | 0:813b917360ac | 548 | |
| higedura | 0:813b917360ac | 549 | for( int i=0;i<N;i++){ yBuf[i] = y[i]; } |
| higedura | 0:813b917360ac | 550 | pk1 = func (yBuf,x); |
| higedura | 0:813b917360ac | 551 | for( int i=0;i<N;i++ ){ k[i][0] = *pk1; pk1 = pk1+1; } |
| higedura | 0:813b917360ac | 552 | |
| higedura | 0:813b917360ac | 553 | for( int i=0;i<N;i++){ yBuf[i] = y[i]+0.5*dt*k[i][1]; } |
| higedura | 0:813b917360ac | 554 | pk2 = func (yBuf,x); |
| higedura | 0:813b917360ac | 555 | for( int i=0;i<N;i++ ){ k[i][1] = *pk2; pk2 = pk2+1; } |
| higedura | 0:813b917360ac | 556 | |
| higedura | 0:813b917360ac | 557 | for( int i=0;i<N;i++){ yBuf[i] = y[i]+0.5*dt*k[i][2]; } |
| higedura | 0:813b917360ac | 558 | pk3 = func (yBuf,x); |
| higedura | 0:813b917360ac | 559 | for( int i=0;i<N;i++ ){ k[i][2] = *pk3; pk3 = pk3+1; } |
| higedura | 0:813b917360ac | 560 | |
| higedura | 0:813b917360ac | 561 | for( int i=0;i<N;i++){ yBuf[i] = y[i]+dt*k[i][3]; } |
| higedura | 0:813b917360ac | 562 | pk4 = func (yBuf,x); |
| higedura | 0:813b917360ac | 563 | for( int i=0;i<N;i++ ){ k[i][3] = *pk4; pk4 = pk4+1; } |
| higedura | 0:813b917360ac | 564 | |
| higedura | 0:813b917360ac | 565 | for( int i=0;i<N;i++){ y[i] = y[i]+dt*(k[i][0]+2.0*k[i][1]+2.0*k[i][2]+k[i][3])/6.0; } |
| higedura | 0:813b917360ac | 566 | |
| higedura | 0:813b917360ac | 567 | return p_y; |
| higedura | 0:813b917360ac | 568 | |
| higedura | 0:813b917360ac | 569 | } |
| higedura | 0:813b917360ac | 570 | |
| higedura | 0:813b917360ac | 571 | double* func( double y[N], double x[N] ){ |
| higedura | 0:813b917360ac | 572 | |
| higedura | 0:813b917360ac | 573 | double f[N] = {0}; |
| higedura | 0:813b917360ac | 574 | double* p_f = f; |
| higedura | 0:813b917360ac | 575 | |
| higedura | 0:813b917360ac | 576 | f[0] = x[0]+x[1]*sin(y[0])*tan(y[1])+x[2]*cos(y[0])*tan(y[1]); |
| higedura | 0:813b917360ac | 577 | f[1] = x[1]*cos(y[0])-x[2]*sin(y[0]); |
| higedura | 0:813b917360ac | 578 | f[2] = x[1]*sin(y[0])/cos(y[1])+x[2]*cos(y[0])/cos(y[1]); |
| higedura | 0:813b917360ac | 579 | |
| higedura | 0:813b917360ac | 580 | return p_f; |
| higedura | 0:813b917360ac | 581 | |
| higedura | 0:813b917360ac | 582 | } |
| higedura | 0:813b917360ac | 583 | |
| higedura | 0:813b917360ac | 584 | double* calib(){ |
| higedura | 0:813b917360ac | 585 | |
| higedura | 0:813b917360ac | 586 | int calib_loop = 1000; |
| higedura | 0:813b917360ac | 587 | |
| higedura | 0:813b917360ac | 588 | int bit_acc[N] = {0}; // Buffer of the accelerometer |
| higedura | 0:813b917360ac | 589 | int get_mag[N] = {0}; // Buffer of the compass |
| higedura | 0:813b917360ac | 590 | |
| higedura | 0:813b917360ac | 591 | double calib_acc[N] = {0}; |
| higedura | 0:813b917360ac | 592 | double calib_gyro_buf[N] = {0}; |
| higedura | 0:813b917360ac | 593 | double calib_gyro[N] = {0}; |
| higedura | 0:813b917360ac | 594 | double compass_basis_buf[N] = {0}; |
| higedura | 0:813b917360ac | 595 | double compass_basis_rad = 0; |
| higedura | 0:813b917360ac | 596 | double calib_result[7] = {0}; |
| higedura | 0:813b917360ac | 597 | |
| higedura | 0:813b917360ac | 598 | double* p_calib_result = calib_result; |
| higedura | 0:813b917360ac | 599 | |
| higedura | 0:813b917360ac | 600 | pc.printf("Don't touch... Calibrating now!!\n\r"); |
| higedura | 0:813b917360ac | 601 | xbee.printf("Don't touch... Calibrating now!!\n\r"); |
| higedura | 0:813b917360ac | 602 | led1 = 1; |
| higedura | 0:813b917360ac | 603 | |
| higedura | 0:813b917360ac | 604 | for( int i=0;i<calib_loop;i++ ){ |
| higedura | 0:813b917360ac | 605 | |
| higedura | 0:813b917360ac | 606 | accelerometer.getOutput(bit_acc); |
| higedura | 0:813b917360ac | 607 | compass.readData(get_mag); |
| higedura | 0:813b917360ac | 608 | |
| higedura | 0:813b917360ac | 609 | calib_gyro_buf[0] = gyro.getGyroX(); |
| higedura | 0:813b917360ac | 610 | calib_gyro_buf[1] = gyro.getGyroY(); |
| higedura | 0:813b917360ac | 611 | calib_gyro_buf[2] = gyro.getGyroZ(); |
| higedura | 0:813b917360ac | 612 | |
| higedura | 0:813b917360ac | 613 | for( int j=0;j<N;j++ ){ |
| higedura | 0:813b917360ac | 614 | calib_acc[j] += (int16_t)bit_acc[j]; |
| higedura | 0:813b917360ac | 615 | calib_gyro[j] += calib_gyro_buf[j]; |
| higedura | 0:813b917360ac | 616 | compass_basis_buf[j] += (int16_t)get_mag[j]; |
| higedura | 0:813b917360ac | 617 | } |
| higedura | 0:813b917360ac | 618 | |
| higedura | 0:813b917360ac | 619 | if( i>calib_loop*3/4 ){ |
| higedura | 0:813b917360ac | 620 | led4 = 1; |
| higedura | 0:813b917360ac | 621 | }else if( i>calib_loop/2 ){ |
| higedura | 0:813b917360ac | 622 | led3 = 1; |
| higedura | 0:813b917360ac | 623 | }else if( i>calib_loop/4 ){ |
| higedura | 0:813b917360ac | 624 | led2 = 1; |
| higedura | 0:813b917360ac | 625 | } |
| higedura | 0:813b917360ac | 626 | |
| higedura | 0:813b917360ac | 627 | } |
| higedura | 0:813b917360ac | 628 | |
| higedura | 0:813b917360ac | 629 | for( int i=0;i<N;i++ ){ |
| higedura | 0:813b917360ac | 630 | calib_acc[i] = calib_acc[i]/calib_loop; |
| higedura | 0:813b917360ac | 631 | calib_gyro[i] = calib_gyro[i]/calib_loop; |
| higedura | 0:813b917360ac | 632 | compass_basis_buf[i] = compass_basis_buf[i]/calib_loop; |
| higedura | 0:813b917360ac | 633 | } |
| higedura | 0:813b917360ac | 634 | |
| higedura | 0:813b917360ac | 635 | compass_basis_rad = compass_basis_buf[1]/compass_basis_buf[0]; |
| higedura | 0:813b917360ac | 636 | compass_basis_rad = atan(compass_basis_rad); |
| higedura | 0:813b917360ac | 637 | led1 = 0; led2 = 0; led3 = 0; led4 = 0; |
| higedura | 0:813b917360ac | 638 | |
| higedura | 0:813b917360ac | 639 | pc.printf(" accX accY accZ gyroX gyroY gyroZ yaw_basis\n\r"); |
| higedura | 0:813b917360ac | 640 | xbee.printf(" accX accY accZ gyroX gyroY gyroZ yaw_basis\n\r"); |
| higedura | 0:813b917360ac | 641 | pc.printf("%6.1f, %6.1f, %6.1f, %6.1f, %6.1f, %6.1f, %6.1f\n\r",calib_acc[0],calib_acc[1],calib_acc[2],calib_gyro[0],calib_gyro[1],calib_gyro[2],compass_basis_rad*180/pi); |
| higedura | 0:813b917360ac | 642 | xbee.printf("%6.1f, %6.1f, %6.1f, %6.1f, %6.1f, %6.1f, %6.1f\n\r",calib_acc[0],calib_acc[1],calib_acc[2],calib_gyro[0],calib_gyro[1],calib_gyro[2],compass_basis_rad*180/pi); |
| higedura | 0:813b917360ac | 643 | |
| higedura | 0:813b917360ac | 644 | for( int i=0;i<3;i++ ){ calib_result[i] = calib_acc[i]; } |
| higedura | 0:813b917360ac | 645 | for( int i=3;i<6;i++ ){ calib_result[i] = calib_gyro[i-3]; } |
| higedura | 0:813b917360ac | 646 | calib_result[6] = compass_basis_rad; |
| higedura | 0:813b917360ac | 647 | |
| higedura | 0:813b917360ac | 648 | if( calib_result[0]==0 && calib_result[1]==0 && calib_result[2]==0 ){ |
| higedura | 0:813b917360ac | 649 | pc.printf("Accelerometer is not available.\r\n"); |
| higedura | 0:813b917360ac | 650 | xbee.printf("Accelerometer is not available.\r\n"); |
| higedura | 0:813b917360ac | 651 | } |
| higedura | 0:813b917360ac | 652 | |
| higedura | 0:813b917360ac | 653 | for( int i=0;i<3;i++ ){ |
| higedura | 0:813b917360ac | 654 | wait(.5); |
| higedura | 0:813b917360ac | 655 | led1 = 1; led2 = 1; led3 = 1; led4 = 1; |
| higedura | 0:813b917360ac | 656 | wait(.5); |
| higedura | 0:813b917360ac | 657 | led1 = 0; led2 = 0; led3 = 0; led4 = 0; |
| higedura | 0:813b917360ac | 658 | } |
| higedura | 0:813b917360ac | 659 | |
| higedura | 0:813b917360ac | 660 | return p_calib_result; |
| higedura | 0:813b917360ac | 661 | |
| higedura | 0:813b917360ac | 662 | } |