Полный текст программы реализации модулей статистической обработки анализатора "Тензотрем"

// diplom. cpp: определяет точку входа для консольного приложения.

#include "stdafx. h"

#include <stdio. h>

#include <malloc. h>

// Транспонирование матрицы

void transponse (double ** x, double ** xt)

{

for (int i = 0; i<512; i++)

for (int j = 0; j<512; j++)

xt [i] [j] = x [j] [i];

}

// структура complex для работы с комплексными числами

// re - действительная часть числа

// im - мнимая часть числа

typedef struct complex

{

complex (double real, double imag)

{

re = real;

im = imag;

}

double re;

doubleim;

};

// Функции работы с комплексными числами

// умножение

complex mul (complex a, complex b)

{return complex (a. re*b. re - a. im*b. im, a. re*b. im + a. im*b. re); }

// сложение

complex add (complex a, complex b)

{return complex (a. re + b. re, a. im + b. im); }

// вычитание

complex sub (complex a, complex b)

{return complex (a. re - b. re, a. im - b. im); }

// вычисление степени числа

double powx (double x, int k)

{

double res = 1;

for (int i = 0; i < k; i++)

res = res*x;

return res;

}

// вычисление факториала

intfact (intn)

{

int res = 1;

for (int i = 1; i < n; i++)

res = res*i;

return res;

}

// разложение синуса вряд тейлора

doublesinx (doublex)

{

double res = 0;

for (int i = 0; i < 100; i++)

res = powx (-1, i-1) *powx (x, 2*i-1) /fact (2*i-1) + res;

returnres;

}

// БПФ

// Аргументы:

// mas - указатель на массив исходных данных

// n - размерность массива mas (является степенью 2)

// y - указатель на массив выходных данных размерности n

void FFT (double * mas, int n, complex * y)

{

double pi = 3.1415926535;

if (n==1)

{

y = &complex (mas [0], 0);

return;

}

complex wn (sinx (2*pi/n + pi/2), sinx (2*pi/n));

complex w (1, 0);

int i = 0, j = n/2, k = 0;

double *a0 = (double *) malloc (sizeof (double) * (n/2)),

*a1 = (double *) malloc (sizeof (double) * (n/2));

while (k < j-1)

{

a0 [k] = mas [i];

a1 [k] = mas [i+1];

k++;

i = i+2;

}

complex *y0 = (complex *) malloc (sizeof (complex) * (n/2));

complex *y1 = (complex *) malloc (sizeof (complex) * (n/2));

FFT (&a0 [0], j, &y0 [0]);

FFT (&a1 [0], j, &y1 [0]);

for (k = 0; k < j; k++)

{

y [k] = add (y0 [k], mul (w,y1 [k]));

y [k+j] = sub (y0 [k], mul (w,y1 [k]));

w = mul (w, wn);

}

}

// Корреляция

// Аргументы:

// vect_1 - указатель на первый массив исходных данных

// vect_2 - указатель на второй массив исходных данных

// corr - указатель на массив выходных данных

// n - размерность данных

// примечание: предполагается, что размерности vect_1 и vect_2 // совпадают

void correlation (double *vect_1, double *vect_2, double *corr, int n)

{

int i, j;

for (i = 0; i < n; i++)

for (j = 0; j<n - i - 1; j++) // т к в выборке всего n

// значений, а индекс i+j не должен становиться больше

// n, то второй цикл до n - i

corr [i] = corr [i] + vect_1 [j] *vect_2 [i+j] /n;

}

// Автокорреляция

// Аргументы:

// vect - указатель на массив исходных данных

// n - размерность массива mas (явялется степенью 2)

// acorr - указатель на массив выходных данных

void autocorrelation (double * vect, int n, double * acorr)

{

int i, j;

for (i = 0; i < n; i++)

for (j = 0; j<n - i - 1; j++) // т к в выборке всего

// n значений, а индекс i+j не должен становиться

// больше n, то второй цикл до n-i

acorr [i] = acorr [i] + vect [j] *vect [i+j] /n;

}

// умножение матриц размера nxm

// Аргументы:

// А - указатель на первый множитель

// В - указатель на второй множитель

// n, m - размерность матриц

// R - указатель на результат

void mulMatr (double ** A, double ** B, int n, int m, double ** R)

{

int i, j, k;

for (i = 0; i < n; i++)

for (j = 0; j < m; j++)

for (k = 0; k < m; k++)

R [i] [j] = R [i] [j] + A [i] [k] *B [k] [j];

}

// Вычисление корня числа

// Аргументы:

// x - число, квадратный корень которого надо вычислить

// Примечание: вычисление ведется методом Ньютона

doublesqrtx (doublex)

{

int i;

double res = 1;

for (i = 0; i < 100; i++)

res = (res + x/res) /2;

return res;

}

// QR-алгоритм

// Аргументы:

// A - указатель на входную квадратную матрицу nxn

// n - размерность данных

// R - указатель на результирующую матрицу с СЧ матрицы А на

// главной диагонали размерности nxn

// Q - указатель на результирующую матрицу СВ матрицы А

// размерности nxn

void QR (double ** A, int n, double ** R, double ** Q)

{

double **E = (double **) malloc (sizeof (double *) * n);

double **T = (double **) malloc (sizeof (double *) * n);

int i, j, k, g;

for (i = 0; i < n; i++)

{

E [i] = (double *) malloc (sizeof (double) * n);

T [i] = (double *) malloc (sizeof (double) * n);

}

for (i = 0; i < n; i++)

{

for (j = 0; j < n; j++)

{

R [i] [j] = A [i] [j];

Q [i] [j] = 1;

if (i == j)

E [i] [j] = 1;

else

E [i] [j] = 0;

}

}

for (k = 0; k < (n-1); k++)

{

for (i = k+1; i < n; i++)

{

for (j = 0; j < n; j++)

{

for (g = 0; g < n; g++)

{

T [j] [g] = E [j] [g];

}

}

if (R [i] [k]! = 0)

{

T [k] [k] = R [k] [k] /sqrtx (powx (R [k] [k],2) +powx (A [i] [k],2));

T [i] [i] = T [k] [k];

T [k] [i] = R [i] [k] /sqrtx (powx (R [k] [k],2) +powx (A [i] [k],2));

T [i] [k] = T [k] [i];

}

mulMatr (&T [0], &R [0], n, n, &R [0]);

mulMatr (&T [0], &Q [0], n, n, &Q [0]);

}

}

free (E);

free (T);

}

// Трансспонирование матрицы размера nxm

// Аргументы:

// x - указатель на исходную матрицу

// xt - указатель на результирующую матрицу

// n,m - размерность матриц

void transp (double **x, double **xt, int n, int m)

{

int i, j;

for (i = 0; i < n; i++)

for (j = 0; j < m; j++)

xt [i] [j] = x [j] [i];

}

// Сортировка СЧ и СВ

// Аргументы:

// eig - указатель на массив СЧ

// Q - указатель на массив СВ

// n - размерность массивов

void sort (double *eig, double ** Q, int n)

{

int i, j;

bool fl = true;

double a;

while (fl)

{

fl = false;

for (i = 0; i < n-1; i++)

if (eig [i] < eig [i+1])

{

fl = true;

a = eig [i];

eig [i] = eig [i+1];

eig [i+1] = a;

for (j = 0; j < n-1; j++)

{

a = Q [i] [j];

Q [i] [j] = Q [i+1] [j];

Q [i+1] [j] = a;

}

}

}

}

// Метод главных компонент

// Аргументы:

// vect - указатель на одномерный массив исходных данных

// y - указатель на одномерный массив результата

// n - размерность входных данных

// k - размерность результата

void PCA (double *vect, double *y, int n, int k)

{

inti, j, m = n/k, h;

// так как vect - вектор-строка, а алгоритм расчитан на

// вектор=столбец, сформируем вектор х=transp (vect)

double **x = (double **) malloc (sizeof (double*));

x [0] = (double *) malloc (sizeof (double) * n);

for (i = 0; i < n; i++)

x [0] [i] = vect [i];

double **Xm = (double **) malloc (sizeof (double *) * n);

double **XT = (double **) malloc (sizeof (double *) * m);

for (i = 0; i < n; i++)

Xm [i] = (double *) malloc (sizeof (double) * m);

for (i = 0; i < m; i++)

XT [i] = (double *) malloc (sizeof (double) * n);

double **Rxx = (double **) malloc (sizeof (double *) * n);

for (i = 0; i < n; i++)

Rxx [i] = (double *) malloc (sizeof (double) * n);

// сформируем матрицу Х исходных векторов

// размерность матрицы определим как 1000хn/1000

h = 0;

for (i = 0; i < m; i++)

for (j = 0; j < n; j++)

{

Xm [j] [i] = x [0] [h];

h++;

}

transp (&Xm [0], &XT [0], n, m);

mulMatr (&Xm [0], &XT [0], n, m, &Rxx [0]);

for (i = 0; i < n; i++)

for (j = 0; j < n; j++)

Rxx [i] [j] = Rxx [i] [j] /m;

double **R = (double **) malloc (sizeof (double *) * n);

double **Q = (double **) malloc (sizeof (double *) * n);

for (i = 0; i < n; i++)

{

R [i] = (double *) malloc (sizeof (double) * n);

Q [i] = (double *) malloc (sizeof (double) * n);

}

QR (&Rxx [0], n, &R [0], &Q [0]);

double *eig = (double *) malloc (sizeof (double) * n);

for (i = 0; i < n; i++)

eig [i] = R [i] [i];

sort (&eig [0], &Q [0], n);

double **W = (double **) malloc (sizeof (double*) * k);

for (i = 0; i < k; i++)

W [i] = (double *) malloc (sizeof (double) * n);

for (i = 0; i < k; i++)

for (j = 0; j < n; j++)

W [i] [j] = Q [j] [i];

for (i = 0; i < k; i++)

for (j = 0; j < n; j++)

y [i] = W [i] [j] *x [0] [j];

free (x);

free (Xm);

free (XT);

free (Rxx);

free (W);

}

// Построение графиков

// Аргументы:

// vect - указатель на массив значений

// n - количество значений

#include<FL/FL. H>

#include "Simple_window. h"

#include "Graph. h"

void plot (double * vect, int n)

{

using namespace Graph_lib;

const int xmax = 600; // размер окна

const int ymax = 400;

const int x_orig = xmax/2; // точка (0,0) - центр окна

const int y_orig = ymax/2;

const Point orig (x_orig, y_orig);

const int r_min = - 1000; // диапазон

const int r_max = 1000;

const int n_noits = (r_max - r_min) *n; // количество точек в

// диапазоне

constintx_scale = 30; // масштабирующий множитель

const int y_scale = 30;

Point p1 (100,100);

Simple_window win (p1, 600, 400, "Plot"); // Создание окна

const int xlen = xmax - 40; // создаем оси

constintylen = ymax - 40; // делаем их чуть меньше

// размера окна

Axisx (Axis:: x, Point (20, y_orig), xlen, xlen/x_scale, "x");

// создание координатных осей

Axisy (Axis:: y, Point (x_orig, ylen + 20), ylen, ylen/y_scale, "y");

// создание координатных осей

win. attach (x);

win. attach (y);

double max = vect [0];

int i;

for (i = 0; i < n; i++)

if (vect [i] > max)

max = vect [i];

int * y1 = (int*) malloc (sizeof (int) * n);

for (i = 0; i < n; i++)

y1 [i] = vect [i] *ymax/max;

Open_polyline poly;

double t = r_min;

double dt = (r_max - r_min) /n;

for (i = 0; i < n; i++)

{

poly. add (Point (t, y1 [i]));

t = t+dt;

}

poly. set_color (Color:: black);

win. attach (poly); // добавление poly в окно

win. wait_for_button (); // вывод окна

return;

}

// Чтение из файла

// Аргументы:

// vect - указатель на вектор значений, считываемых из файла

#include <stdlib. h>

int readFile (double * vect)

{

int i;

printf ("%s","filename: ");

char filename [200]; // ={"I: \\1. txt"};

scanf ("%s", &filename);

FILE * file;

file = fopen (filename, "r");

i = 0;

if (file! = NULL)

{

while (true)

{

vect = (double *) realloc (vect, sizeof (double) * (i+1));

if (vect! = NULL)

{

vect [i] = fscanf (file,"%f",&vect [i]); // читаем из файла

// очередное число

}

if (vect [i] == EOF || vect == NULL)

{

fclose (file);

return i;

}

i++;

}

return i;

}

}

// Главная функция

int _tmain (int argc, _TCHAR* argv [])

{

bool f0 = true, f1 = true;

while (f0)

{

f1 = true;

// ввод имени файла и данных из него

printf ("%s", "Please enter the file name\n");

double * vect = (double *) malloc (sizeof (double));

int n = readFile (vect);

// выбор метода обработки

int num;

printf ("%s", "Please select the method\n");

printf ("%s", "1. FFT\n");

printf ("%s", "2. PCA\n");

printf ("%s", "3. Correlation\n");

printf ("%s", "4. Autocorrelation\n");

printf ("%s", "5. To attach another file\n");

printf ("%s", "6. Exit\n");

while (f1)

{

printf ("%s", "Number: ");

scanf ("%d", &num);

if (num == 1)

{

inti = n;

if (n& (n-1)) // если побитовое и дает 1, то n // не // равно степени 2

{

while (i& (i-1))

{

vect = (double *) realloc (vect, sizeof (double) * (i+1));

vect [i] = 0;

}

}

complex * y = (complex*) malloc (sizeof (complex) * i);

FFT (vect, i, y);

double * y1 = (double*) malloc (sizeof (double) * i);

for (int j = 0; j < i; j++)

y1 [j] = sqrt (powx (y [j]. re,2) + powx (y [j]. im,2));

plot (y1, n);

free (y);

free (y1);

}

else if (num == 2)

{

double * rPCA = (double*) malloc (sizeof (double));

PCA (vect, rPCA, n, 1000);

plot (rPCA, 1000);

free (rPCA);

}

else if (num == 3)

{

printf ("%s","Please enter the name of the second file\n");

double * vect_1 = (double *) malloc (sizeof (double));

readFile (vect_1);

double * corr = (double *) malloc (sizeof (double) * n);

correlation (vect, vect_1, corr, n);

plot (corr, n);

free (vect_1);

free (corr);

}

else if (num == 4)

{

double * acorr = (double *) malloc (sizeof (double) *n);

autocorrelation (&vect [0], n, &acorr [0]);

plot (acorr, n);

free (acorr);

}

else if (num == 5)

f1 = false;

else if (num == 6)

{

f1 = false;

f0 = false;

}

}

}

return 0;

}

Размещено на Allbest.ru