Example of a least square caculation, Matrix calculation based on CMSIS-DSP library, x(t)=a+b*t+c*t^2.

Dependencies:   dsp mbed

Comment out line 226 of "arm_math.h" to prevent a lot of compiler warnings

main.cpp

Committer:
GerritPathuis
Date:
2015-03-10
Revision:
0:9087ea0ad689

File content as of revision 0:9087ea0ad689:

///////////////////////////////////////
// FRDM-K22F                         //
// Least square calculation (Matrix) //
// Based on the DSP libtary          //
// Digital Signal Processing         //
//                                   //
// Results can be read via Tera Term //
///////////////////////////////////////
#include "mbed.h"
#include "arm_math.h"

#define NUMSAMPLES 51
#define NUMUNKNOWS 3

Serial pc(USBTX, USBRX); // tx, rx

// time is evenly spaced bit this in not required
float32_t tData[NUMSAMPLES] = {
    0.0f,    0.1f,    0.2f,    0.3f,    0.4f,    0.5f,    0.6f,    0.7f,
    0.8f,    0.9f,    1.0f,    1.1f,    1.2f,    1.3f,    1.4f,    1.5f,
    1.6f,    1.7f,    1.8f,    1.9f,    2.0f,    2.1f,    2.2f,    2.3f,
    2.4f,    2.5f,    2.6f,    2.7f,    2.8f,    2.9f,    3.0f,    3.1f,
    3.2f,    3.3f,    3.4f,    3.5f,    3.6f,    3.7f,    3.8f,    3.9f,
    4.0f,    4.1f,    4.2f,    4.3f,    4.4f,    4.5f,    4.6f,    4.7f,
    4.8f,    4.9f,    5.0f
};

float32_t xData[NUMSAMPLES] = {
    07.4213f,   21.7231f,    -7.2828f,    21.2254f,    20.2221f,    10.3585f,    20.2033f,   29.2690f,
    57.7152f,   53.6075f,    22.8209f,    59.8714f,    43.1712f,    38.4436f,    46.0499f,   39.8803f,
    41.5188f,   55.2256f,    55.1803f,    55.6495f,    49.8920f,    34.8721f,    50.0859f,   57.0099f,
    47.3032f,   50.8975f,    47.4671f,    38.0605f,    41.4790f,    31.2737f,    42.9272f,   24.6954f,
    23.1770f,   22.9120f,    3.2977f,     35.6270f,    23.7935f,    12.0286f,    25.7104f,   -2.4601f,
    06.7021f,   01.6804f,    02.0671f,    -2.2891f,    -16.2070f,   -14.2204f,   -20.1780f,  -18.9303f,
    -20.4859f,  -25.8338f,   -47.2892f
};



float32_t AData[NUMSAMPLES * NUMUNKNOWS];
float32_t ATData[NUMSAMPLES * NUMUNKNOWS];
float32_t ATAData[NUMUNKNOWS * NUMUNKNOWS];
float32_t invATAData[NUMUNKNOWS * NUMUNKNOWS];
float32_t BData[NUMSAMPLES * NUMSAMPLES];
float32_t cData[NUMUNKNOWS];

arm_matrix_instance_f32 t=  {NUMSAMPLES, 1, tData};
arm_matrix_instance_f32 x=  {NUMSAMPLES, 1, xData};
arm_matrix_instance_f32 A=  {NUMSAMPLES, NUMUNKNOWS, AData};
arm_matrix_instance_f32 AT= {NUMUNKNOWS, NUMSAMPLES,  ATData};
arm_matrix_instance_f32 ATA= {NUMUNKNOWS, NUMUNKNOWS, ATAData};
arm_matrix_instance_f32 invATA= {NUMUNKNOWS, NUMUNKNOWS, invATAData};
arm_matrix_instance_f32 B=  {NUMUNKNOWS, NUMSAMPLES, BData};
arm_matrix_instance_f32 c=  {NUMUNKNOWS, 1, cData};

main(void)
{
    int i, time;
    float y;

    pc.printf("/n/nMatrix calculation (Least squares) with 3 Unknowns \n\r");
    pc.printf("We calculates a, b and c from; x(t)= a + b*t + c*t*t \n\r");
    pc.printf("Results approx a= 8.701225, b= 38.880058, c= -9.792965\n\r");

    y= sqrtf(xData[0]);
    cData[0]= y;

    for(i=0; i<NUMSAMPLES; i++) {
        AData[i*NUMUNKNOWS +0] = 1.0f;
        AData[i*NUMUNKNOWS +1] = tData[i];
        AData[i*NUMUNKNOWS +2] = tData[i] * tData[i];
    }

    /* calculation of A transpose */
    arm_mat_trans_f32(&A, &AT);
    arm_mat_mult_f32(&AT, &A, &ATA);
    arm_mat_inverse_f32(&ATA, &invATA);
    arm_mat_mult_f32(&invATA, &AT, &B);
    arm_mat_mult_f32(&B, &x, &c);


    pc.printf("Results a= %f, b= %f, c= %f \n\r", cData[0], cData[1], cData[2]);

    for(time=0; time<6; time++) {
        y= cData[0] + cData[1]* time + cData[2] * time * time;
        pc.printf("time= %d, x= %f \n\r", time, y );
    }
}