arm_matrix_example_f32.c
1 /* ---------------------------------------------------------------------- 2 * Copyright (C) 2010-2012 ARM Limited. All rights reserved. 3 * 4 * $Date: 17. January 2013 5 * $Revision: V1.4.0 6 * 7 * Project: CMSIS DSP Library 8 * Title: arm_matrix_example_f32.c 9 * 10 * Description: Example code demonstrating least square fit to data 11 * using matrix functions 12 * 13 * Target Processor: Cortex-M4/Cortex-M3 14 * 15 * Redistribution and use in source and binary forms, with or without 16 * modification, are permitted provided that the following conditions 17 * are met: 18 * - Redistributions of source code must retain the above copyright 19 * notice, this list of conditions and the following disclaimer. 20 * - Redistributions in binary form must reproduce the above copyright 21 * notice, this list of conditions and the following disclaimer in 22 * the documentation and/or other materials provided with the 23 * distribution. 24 * - Neither the name of ARM LIMITED nor the names of its contributors 25 * may be used to endorse or promote products derived from this 26 * software without specific prior written permission. 27 * 28 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS 29 * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT 30 * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS 31 * FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE 32 * COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, 33 * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, 34 * BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; 35 * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER 36 * CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT 37 * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN 38 * ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE 39 * POSSIBILITY OF SUCH DAMAGE. 40 * -------------------------------------------------------------------- */ 41 42 /** 43 * @ingroup groupExamples 44 */ 45 46 /** 47 * @defgroup MatrixExample Matrix Example 48 * 49 * \par Description: 50 * \par 51 * Demonstrates the use of Matrix Transpose, Matrix Muliplication, and Matrix Inverse 52 * functions to apply least squares fitting to input data. Least squares fitting is 53 * the procedure for finding the best-fitting curve that minimizes the sum of the 54 * squares of the offsets (least square error) from a given set of data. 55 * 56 * \par Algorithm: 57 * \par 58 * The linear combination of parameters considered is as follows: 59 * \par 60 * <code>A * X = B</code>, where \c X is the unknown value and can be estimated 61 * from \c A & \c B. 62 * \par 63 * The least squares estimate \c X is given by the following equation: 64 * \par 65 * <code>X = Inverse(A<sup>T</sup> * A) * A<sup>T</sup> * B</code> 66 * 67 * \par Block Diagram: 68 * \par 69 * 70 * \par Variables Description: 71 * \par 72 * \li \c A_f32 input matrix in the linear combination equation 73 * \li \c B_f32 output matrix in the linear combination equation 74 * \li \c X_f32 unknown matrix estimated using \c A_f32 & \c B_f32 matrices 75 * 76 * \par CMSIS DSP Software Library Functions Used: 77 * \par 78 * - arm_mat_init_f32() 79 * - arm_mat_trans_f32() 80 * - arm_mat_mult_f32() 81 * - arm_mat_inverse_f32() 82 * 83 * <b> Refer </b> 84 * \link arm_matrix_example_f32.c \endlink 85 * 86 */ 87 88 89 /** \example arm_matrix_example_f32.c 90 */ 91 92 #include "arm_math.h" 93 #include "math_helper.h" 94 95 #if defined(SEMIHOSTING) 96 #include <stdio.h> 97 #endif 98 99 #define SNR_THRESHOLD 90 100 101 /* -------------------------------------------------------------------------------- 102 * Test input data(Cycles) taken from FIR Q15 module for differant cases of blockSize 103 * and tapSize 104 * --------------------------------------------------------------------------------- */ 105 106 const float32_t B_f32[4] = 107 { 108 782.0, 7577.0, 470.0, 4505.0 109 }; 110 111 /* -------------------------------------------------------------------------------- 112 * Formula to fit is C1 + C2 * numTaps + C3 * blockSize + C4 * numTaps * blockSize 113 * -------------------------------------------------------------------------------- */ 114 115 const float32_t A_f32[16] = 116 { 117 /* Const, numTaps, blockSize, numTaps*blockSize */ 118 1.0, 32.0, 4.0, 128.0, 119 1.0, 32.0, 64.0, 2048.0, 120 1.0, 16.0, 4.0, 64.0, 121 1.0, 16.0, 64.0, 1024.0, 122 }; 123 124 125 /* ---------------------------------------------------------------------- 126 * Temporary buffers for storing intermediate values 127 * ------------------------------------------------------------------- */ 128 /* Transpose of A Buffer */ 129 float32_t AT_f32[16]; 130 /* (Transpose of A * A) Buffer */ 131 float32_t ATMA_f32[16]; 132 /* Inverse(Transpose of A * A) Buffer */ 133 float32_t ATMAI_f32[16]; 134 /* Test Output Buffer */ 135 float32_t X_f32[4]; 136 137 /* ---------------------------------------------------------------------- 138 * Reference ouput buffer C1, C2, C3 and C4 taken from MATLAB 139 * ------------------------------------------------------------------- */ 140 const float32_t xRef_f32[4] = {73.0, 8.0, 21.25, 2.875}; 141 142 float32_t snr; 143 144 145 /* ---------------------------------------------------------------------- 146 * Max magnitude FFT Bin test 147 * ------------------------------------------------------------------- */ 148 149 int32_t main(void) 150 { 151 152 arm_matrix_instance_f32 A; /* Matrix A Instance */ 153 arm_matrix_instance_f32 AT; /* Matrix AT(A transpose) instance */ 154 arm_matrix_instance_f32 ATMA; /* Matrix ATMA( AT multiply with A) instance */ 155 arm_matrix_instance_f32 ATMAI; /* Matrix ATMAI(Inverse of ATMA) instance */ 156 arm_matrix_instance_f32 B; /* Matrix B instance */ 157 arm_matrix_instance_f32 X; /* Matrix X(Unknown Matrix) instance */ 158 159 uint32_t srcRows, srcColumns; /* Temporary variables */ 160 arm_status status; 161 162 /* Initialise A Matrix Instance with numRows, numCols and data array(A_f32) */ 163 srcRows = 4; 164 srcColumns = 4; 165 arm_mat_init_f32(&A, srcRows, srcColumns, (float32_t *)A_f32); 166 167 /* Initialise Matrix Instance AT with numRows, numCols and data array(AT_f32) */ 168 srcRows = 4; 169 srcColumns = 4; 170 arm_mat_init_f32(&AT, srcRows, srcColumns, AT_f32); 171 172 /* calculation of A transpose */ 173 status = arm_mat_trans_f32(&A, &AT); 174 175 176 /* Initialise ATMA Matrix Instance with numRows, numCols and data array(ATMA_f32) */ 177 srcRows = 4; 178 srcColumns = 4; 179 arm_mat_init_f32(&ATMA, srcRows, srcColumns, ATMA_f32); 180 181 /* calculation of AT Multiply with A */ 182 status = arm_mat_mult_f32(&AT, &A, &ATMA); 183 184 /* Initialise ATMAI Matrix Instance with numRows, numCols and data array(ATMAI_f32) */ 185 srcRows = 4; 186 srcColumns = 4; 187 arm_mat_init_f32(&ATMAI, srcRows, srcColumns, ATMAI_f32); 188 189 /* calculation of Inverse((Transpose(A) * A) */ 190 status = arm_mat_inverse_f32(&ATMA, &ATMAI); 191 192 /* calculation of (Inverse((Transpose(A) * A)) * Transpose(A)) */ 193 status = arm_mat_mult_f32(&ATMAI, &AT, &ATMA); 194 195 /* Initialise B Matrix Instance with numRows, numCols and data array(B_f32) */ 196 srcRows = 4; 197 srcColumns = 1; 198 arm_mat_init_f32(&B, srcRows, srcColumns, (float32_t *)B_f32); 199 200 /* Initialise X Matrix Instance with numRows, numCols and data array(X_f32) */ 201 srcRows = 4; 202 srcColumns = 1; 203 arm_mat_init_f32(&X, srcRows, srcColumns, X_f32); 204 205 /* calculation ((Inverse((Transpose(A) * A)) * Transpose(A)) * B) */ 206 status = arm_mat_mult_f32(&ATMA, &B, &X); 207 208 /* Comparison of reference with test output */ 209 snr = arm_snr_f32((float32_t *)xRef_f32, X_f32, 4); 210 211 /*------------------------------------------------------------------------------ 212 * Initialise status depending on SNR calculations 213 *------------------------------------------------------------------------------*/ 214 status = (snr < SNR_THRESHOLD) ? ARM_MATH_TEST_FAILURE : ARM_MATH_SUCCESS; 215 216 if (status != ARM_MATH_SUCCESS) 217 { 218 #if defined (SEMIHOSTING) 219 printf("FAILURE\n"); 220 #else 221 while (1); /* main function does not return */ 222 #endif 223 } 224 else 225 { 226 #if defined (SEMIHOSTING) 227 printf("SUCCESS\n"); 228 #else 229 while (1); /* main function does not return */ 230 #endif 231 } 232 233 } 234 235 /** \endlink */