Svenska ElektronikForumet

Kod: Markera allt

matrix* cut(matrix* a, int start_n, int stop_n, int start_m, int stop_m) {

	int n = a->row;
	float* data = a->data;

	// Create the output
	matrix* out = initMatrix(stop_n - start_n + 1, stop_m - start_m + 1);
	float* ptr = out->data;

	for (int i = start_n; i < stop_n + 1; i++) {
		memcpy(ptr++, data + i*n + start_m, sizeof(stop_m+1));
	}

	return out;
}

Kod: Markera allt

matrix* mul(matrix* a, matrix* b, bool elementWise) {
	// Get the dimensions - they both are the same
	int row_a = a->row;
	int column_a = a->column;
	//int row_b = b->row;
	int column_b = b->column;

	// Data
	float* data_a = a->data;
	float* data_b = b->data;

	if (elementWise == true) {
		matrix* out = initMatrix(row_a, column_a); // Same dimensions for a and b matrix
		float* ptr = out->data;
		// Dot multiply all values
		if (column_b > 1) { // If matrix b is a matrix
			for (int i = 0; i < row_a; i++) {
				for (int j = 0; j < column_a; j++) {
					// Do element wise mutiplication. In MATLAB it is A.*A
					*(ptr++) = *(data_a++) * *(data_b++);
				}
			}
		} else {
			// If matrix b is a vector
			for (int i = 0; i < row_a; i++) {
				for (int j = 0; j < column_a; j++) {
					// Do element wise mutiplication. In MATLAB it is A.*b
					*(ptr++) = *(data_a++) * *(data_b + i);
				}
			}
		}

		return out;
	} else {
		// Do regular mutiplication. In MATLAB it is A*A
		matrix* out = initMatrix(row_a, column_b);
		float* ptr = out->data;

		// Let's take our a matrix
		for (int i = 0; i < a->row; i++) {

			// Then we go through every column of b
			for (int j = 0; j < b->column; j++) {
				data_a = &a->data[i * a->column];
				data_b = &b->data[j];

				*ptr = 0; // Reset
				// And we multiply rows from a with columns of b
				for (int k = 0; k < a->column; k++) {
					*ptr += *data_a * *data_b;
					data_a++;
					data_b += b->column;
				}
				ptr++;
			}
		}

		return out;
	}

}

Kod: Markera allt

>> tic; [u,s,v]=svd(e); toc
Elapsed time is 0.00116205 seconds.
>>

Kod: Markera allt

#include <stdio.h>
#include <stdlib.h>
#include "LinearAlgebra/declareFunctions.h"
#include <time.h>

int main() {

	/*
	 * G(s) = 3.2/(2s^2 + 0.7s + 3.1) - Model
	 */

	float input[72] = {  5,   5,   5,   5,   5,   5,   5,   5,   5,   5,   5,   5,   5,   5,   5,   5,   5,   5,   5,   5,   5,   5,   5,   5,   5,   5,   5,   5,

			   5,   5,   5,   5,   5,   5,   5,   5,   5,   5,   5,   5,   5,   5,   5,   5,   5,   5,   5,   5,   5,   5,   5,   5,   5,   5,   5,   5,

			   5,   5,   5,   5,   5,   5,   5,   5,   5,   5,   5,   5,   5,   5,   5,   5};


	float output[72] = { 0.00000,   1.89649,   3.30777,   3.40331,   3.23615,   3.18151,   3.19049,   3.20013,   3.20132,   3.20031,   3.19990,


			   3.19993,   3.20000,   3.20001,   3.20000,   3.20000,   3.20000,   3.20000,   3.20000,   3.20000,   3.20000,   3.20000,


			   3.20000,   3.20000,   3.20000,   3.20000,   3.20000,   3.20000,   3.20000,   3.20000,   3.20000,   3.20000,   3.20000,


			   3.20000,   3.20000,   3.20000,   3.20000,   3.20000,   3.20000,   3.20000,   3.20000,   3.20000,   3.20000,   3.20000,


			   3.20000,   3.20000,   3.20000,   3.20000,   3.20000,   3.20000,   3.20000,   3.20000,   3.20000,   3.20000,   3.20000,


			   3.20000,   3.20000,   3.20000,   3.20000,   3.20000,   3.20000,   3.20000,   3.20000,   3.20000,   3.20000,   3.20000,


			   3.20000,   3.20000,   3.20000,   3.20000,   3.20000,   3.20000};


	// Create toeplitz matrix
	matrix* toe = toeplitz(input, 72);

	// Create upper triangular matrix of the toeplitz matrix
	matrix* tru = triu(toe, 0);

	// Find the inverse of tru
	matrix* iv = inv(tru);

	// Create vector the horizon - Important! Else, we cannot find the markov-parameters g
    matrix* Y = create(output, 1, 72);

	// Multiply Y with the iv matrix - Find the markov-parameters g
    matrix* G = mul(Y, iv, false);

    // Detete some mats
    freeMatrix(toe);
    freeMatrix(tru);
    freeMatrix(iv);
    freeMatrix(Y);

    /*--------------------------------------*/

    // turn vector G into a normal vector because the function hankel only want 1D vector
    float g[72];
    for(int i = 0; i < 72; i++)
    	g[i] = *(G->data + i);

    // Delete G
    freeMatrix(G);

    // Create hankel matrix H0 and H1
    matrix* H0 = hankel(g, 72, 0);
    matrix* H1 = hankel(g, 72, 1);

    // Cut H1 and H2 to the half - Remember indexing from zero here!
    matrix* H0_half = cut(H0, 0, 35, 0, 35); // 36x36 from 72x72

    // Remove H0
    freeMatrix(H0);

    // Measure the size

    /*
    int n;
    int m;
    size(H0_half, &n, &m);
    printf("The H0_half have the size %dx%d\n\n,", n, m);
    */


    // Do SVD on
    matrix* u =  initMatrix(36,36); // We know the size from the size() command above
    matrix* s =  initMatrix(36,36);
    matrix* v =  initMatrix(36,36);
    clock_t start, end;
    float cpu_time_used;

    start = clock();
    svd(H0_half, u, s, v);
    end = clock();
    cpu_time_used = ((float) (end - start)) / CLOCKS_PER_SEC;
    printf("Total speed of SVD was %f,", cpu_time_used);


    // Delete H1 and H0_half
    freeMatrix(H1);
    freeMatrix(H0_half);
    freeMatrix(u);
    freeMatrix(s);
    freeMatrix(v);


	return EXIT_SUCCESS;
}