39413

Реализация и исследование быстрого алгоритма двумерного вещественного ДПФ по основанию 4 представлением данных в гиперкомплексной алгебре

Курсовая

Информатика, кибернетика и программирование

Заданный алгоритм был реализован программно с помощью технологии Microsoft. NET Framework на языке программирования C++. Написанное приложение состоит из двух сборок: библиотеки классов FFT, содержащей все необходимое для вычисления ДПФ по формуле и БПФ.

Русский

2013-10-04

294.73 KB

20 чел.

Федеральное агентство по образованию

Государственное образовательное учреждение

высшего профессионального образования

«Самарский государственный аэрокосмический

университет  имени академика С.П. Королева»

Кафедра геоинформатики

ПОЯСНИТЕЛЬНАЯ ЗАПИСКА

к курсовому проекту по дисциплине

"Компьютерная алгебра"

Вариант № 11

Cтудент: Терентьева Д.А.
группа 65
7

Руководитель проекта: Чичева М. А.

Самара 2010

Задание

Реализация  и исследование быстрого алгоритма двумерного вещественного ДПФ по основанию 4 представлением данных в гиперкомплексной алгебре.


Содержание

1 Постановка задачи 4

2 Теоретический материал 5

2.1 Двумерное преобразование Фурье 5

2.1.1 Определение двумерного ДПФ в алгебре комплексных чисел  и в четырехмерных алгебрах 5

2.1.2 Алгоритм двумерного ДПФ по основанию 2 6

2.1.3 Алгоритм двумерного ДПФ по основанию 4 7

2.2 Алгебра гиперкомплексных чисел 7

2.3 Прямая сумма комплексных алгебр 9

3 Описание разработанного алгоритма 10

3.1 Описание работы программы 10

4 Исследование алгоритма 12

Приложение А. Текст программы 13

Приложение Б. Принцип работы программы 19


1 Постановка задачи

Требуется разработать и реализовать быстрый алгоритм двумерного преобразования Фурье в гиперкомплексной алгебре по основанию 4.

На входе – квадратная матрица размера , где где .

Тестирование полученной реализации алгоритма, ее исследование и сравнение с алгоритмом двумерного ДПФ.


2 Теоретический материал

2.1 Двумерное преобразование Фурье

2.1.1 Определение двумерного ДПФ в алгебре комплексных чисел
и в четырехмерных алгебрах

Пусть - входной массив размера N×N, тогда его комплексный спектр Фурье

 ,  , (2.1)

где - комплексный корень из единицы степени N.

Под дискретным преобразованием Фурье со специальным представлением данных в рамках курсового проекта будем понимать преобразование вида

 . (2.2)

Основная идея такого преобразования заключается в погружении входного сигнала и корней в четырехмерную алгебру (в рамках курсового проекта это одна из алгебраических структур, описанных в разделе 1). При этом корни лежат в разных экземплярах поля комплексных чисел, вложенных в , у них разные мнимые единицы - i и j.

Отметим, что умножения на степени и в соотношении (2.2) записаны слева и справа от входного сигнала. Это сделано для унификации записи. При использовании алгебры кватернионов это принципиально, так как умножение в ней некоммутативно. При использовании гиперкомплексной алгебры порядок сомножителей значения не имеет.

Поскольку элемент алгебры (кватернион или гиперкомплексное число ) определяется набором четырех вещественных чисел , то комплексный спектр (2.1) может быть получен из спектра (2.2) в алгебре   следующим образом:

 

где - вектор компонент спектра,

 .

Таким образом, мультипликативная сложность вычисления совпадает со сложностью вычисления спектра в алгебре , т.к. умножения на матрицы L, I не требуют выполнения нетривиальных операций вещественного умножения.

Далее излагаются схемы декомпозиции двумерного дискретного преобразования Фурье (2.2) и способы его вычисления.

2.1.2 Алгоритм двумерного ДПФ по основанию 2

Пусть . Представим (2.2) в виде четырех сумм, разделяя входную последовательность по четным и нечетным значениям каждого индекса :

  (2.3)

где

Вычисление спектра для остальных значений пар производится без дополнительных умножений и может быть записано в матричной форме

 

или в следующем виде

 , (2.4)
.


2.1.3 Алгоритм двумерного ДПФ по основанию 4

Разобьем теперь входную последовательность на 16 подпоследовательностей и запишем (2.2) в виде

  (2.5)

где

 , (2.6)

При этом значения спектра для остальных значений аргументов вычисляются без дополнительных умножений, а именно:

 , (2.7)
.

Умножения на степени базовых элементов i и j тривиальны, они сводятся к перестановкам элементов кода и/или смене знака компонент.

2.2 Алгебра гиперкомплексных чисел

Коммутативной гиперкомплексной алгеброй будем называть четырехмерную ассоциативную алгебру над

  (2.2.11)

с определяющими соотношениями для умножений базисных элементов :

  (2.2.12)

Операция сложения в данной алгебре осуществляется покомпонентно. Соотношение (1.12) для умножения базисных элементов индуцирует правило умножения произвольных элементов алгебры

  (2.2.13)

Отметим еще два полезных свойства алгебры :

- поле комплексных чисел  канонически вкладывается в :

 


- справедливы соотношения

 , (2.2.14)
,  .

Непосредственное умножение элементов алгебры по формуле (2.2.13) требует 16 вещественных умножений и 12 вещественных сложений. Как обычно, будем считать, что:

- при оценке вычислительной сложности один из сомножителей предполагается постоянным и, следовательно, все арифметические операции над его компонентами могут быть реализованы заранее;

- умножения на степени числа 2 являются более простыми операциями, чем сложения и умножения, и не учитываются при анализе вычислительной сложности алгоритмов ДПФ.

Пусть - элемент алгебры общего вида; - комплексное число, являющееся константой в контексте рассматриваемых алгоритмов. В соответствии с представлением (2.2.14) вычисление произведения эквивалентно выполнению двух комплексных умножений. При использовании схемы "три умножения, три сложения" искомое произведение принимает вид

  (2.2.15)

и вычисляется посредством шести вещественных умножений и шести вещественных сложений.

Одновременное умножение на два комплексных корня в этой алгебре лучше выполнять по очереди.

Непосредственной проверкой легко убедиться также в справедливости следующего утверждения.

Отображения

  (2.2.16)

сохраняют сумму и произведение элементов алгебры , действуют тождественно на (то есть являются автоморфизмами над ).


2.3 Прямая сумма комплексных алгебр

Известно, что коммутативно-ассоциативная гиперкомплексная алгебра (2.2.11)-(2.2.12) изоморфна прямой сумме комплексных алгебр:

 .

Для перехода от стандартного представления гиперкомплексного числа

  (2.3.1)

к его представлению в прямой сумме комплексных алгебр, выполним замену переменных:

 , , , . (2.3.2)

В этом случае произвольный элемент алгебры примет вид

 

Таблица умножения новых базисных элементов имеет вид, показанный ниже.

Таблица 1 Правила умножения новых базисных элементов

Базисные
элементы

0

0

0

0

0

0

0

0

В таком представлении произведение двух произвольных элементов алгебры

 

разбивается на два независимых произведения, подобных произведениям комплексных чисел

  ,

и требует 6 вещественных умножений и 6 вещественных сложений.

Все свойства гиперкомплексной алгебры в таком представлении сохраняются.


3 Описание разработанного алгоритма

3.1 Описание работы программы

Заданный алгоритм был реализован программно с помощью технологии Microsoft .NET Framework на языке программирования C++.

Написанное приложение состоит из двух сборок: библиотеки классов FFT, содержащей все необходимое для вычисления ДПФ по формуле и БПФ.

Также присутствует консольный модуль, реализующий пользовательский интерфейс:

Рисунок 1 Основное меню программы

После выбора алгоритма пользователь видит консоль с предложением ввести размер квадратного массива:

Рисунок 2 Консоль ввода массива

В результате генерируется файл с подаваемым на вход сигналом, на выходе получаем файл с исходным спектром в комплексном представлении. Далее представлены примеры выполненной работы алгоритмов БПФ и ДПФ для длины преобразования 64:

Рисунок 3 Пример работы алгоритма вычисления БПФ-64

Рисунок 4 Пример работы алгоритма вычисления ДПФ-64


4 Исследование алгоритма

Для наглядного представления приведем таблицу и график сравнения алгоритмов БПФ и ДПФ.

Проведем сравнение результатов непосредственного вычисления дискретного преобразования Фурье и результатов быстрых алгоритмов для длин преобразований 2, 4, 8, 16, 32, 64, 128, 256.

Таблица 2 Сравнение времени выполнения БПФ и  ДПФ

N

Время вычисления БПФ, с

Время вычисления ДПФ, с

8

0,001322

0,001686

16

0,007954

0,018803

32

0,047996

0,455932

64

0,287401

8,510700

128

1,714250

160,6980

256

8,678940

3559,830

 

Рисунок 5 Зависимость времени выполнения алгоритма от длины преобразования

Из Рисунка 5 видно, что быстрый алгоритм Фурье имеет значительное преимущество в скорости перед непосредственным дискретным преобразованием, хотя для небольшого объёма данных различия несущественны. Так же присутствует погрешность, связанная с машинными характеристиками, которая напрямую влияет на скорость алгоритмов.


Приложение А. Текст программы

Class Hypercomplex

#include "StdAfx.h"

#include "Hypercomplex.h"

#include <iostream>

#include <math.h>

#include <float.h>

Hypercomplex::Hypercomplex()

{

 W = 0.0;//1.0

 X = 0.0;

 Y = 0.0;

 Z = 0.0;

}

Hypercomplex::Hypercomplex(const double w, const double x, const double y, const double z)

{

 W = w;

 X = x;

 Y = y;

 Z = z;

}

void Hypercomplex::PrintOn(ostream &strm) const

{

Hypercomplex q(W, X, Y, Z);

 if(abs(q.getW()) < 0.000001)

{

 q.setW(0.0);

}

 if(abs(q.getX()) < 0.000001)

{

 q.setX(0.0);

}

 if(abs(q.getY()) < 0.000001)

{

 q.setY(0.0);

}

 if(abs(q.getZ()) < 0.000001)

{

 q.setZ(0.0);

}

strm << "(" << q.W << "," << q.X << "," << q.Y << "," << q.Z << ")";

}

Hypercomplex operator * (const Hypercomplex &a, const Hypercomplex &b)

{

 double w,x,y,z;

 w = a.W*b.W - (a.X*b.X + a.Y*b.Y + a.Z*b.Z);

 

 x = a.W*b.X + b.W*a.X + a.Y*b.Z + a.Z*b.Y;

 y = a.W*b.Y + b.W*a.Y + a.Z*b.X + a.X*b.Z;

 z = a.W*b.Z + b.W*a.Z + a.X*b.Y + a.Y*b.X;

 return Hypercomplex(w,x,y,z);

}

const Hypercomplex& Hypercomplex::operator *= (const Hypercomplex &q)

{

 double w = W, x = X, y = Y, z = Z;

 double qW = q.W, qX = q.X, qY = q.Y, qZ = q.Z;

W = w*qW - (x*qX + y*qY + z*qZ);

X = w*qX + qW*x + y*qZ + z*qY;

Y = w*qY + qW*y + z*qX + x*qZ;

Z = w*qZ + qW*z + x*qY + y*qX;

 return *this;

}

Hypercomplex operator + (const Hypercomplex &a, const Hypercomplex &b)

{

 double w,x,y,z;

w = a.W+b.W;

x = a.X+b.X;

y = a.Y+b.Y;

z = a.Z+b.Z;

 return Hypercomplex(w,x,y,z);

}

const Hypercomplex& Hypercomplex::operator += (const Hypercomplex &q)

{

W = W+q.W;

X = X+q.X;

Y = Y+q.Y;

Z = Z+q.Z;

 return *this;

}

Hypercomplex operator - (const Hypercomplex &a, const Hypercomplex &b)

{

 double w,x,y,z;

w = a.W-b.W;

x = a.X-b.X;

y = a.Y-b.Y;

z = a.Z-b.Z;

 return Hypercomplex(w,x,y,z);

}

const Hypercomplex& Hypercomplex::operator -= (const Hypercomplex &q)

{

W = W-q.W;

X = X-q.X;

Y = Y-q.Y;

Z = Z-q.Z;

 return *this;

}

Hypercomplex Hypercomplex::operator = (const Hypercomplex &q)

{

W = q.W;

X = q.X;

Y = q.Y;

Z = q.Z;

 return *this;

}

const Hypercomplex& Hypercomplex::operator ~ ()

{

X = -X;

Y = -Y;

Z = -Z;

 return *this;

}

const Hypercomplex& Hypercomplex::exp()

{

 double Mul;

 double Length = sqrt(X*X + Y*Y + Z*Z);

 if (Length > 1.0e-4)

Mul = sin(Length)/Length;

 else

Mul = 1.0;

W = cos(Length);

X *= Mul;

Y *= Mul;

Z *= Mul;

 return *this;

}

const double& Hypercomplex::arg()

{

 return acos(W/sqrt(W*W + X*X + Y*Y + Z*Z));

}

const Hypercomplex& Hypercomplex::sgn()

{

Hypercomplex q = *this;

 return (1/sqrt(W*W + X*X + Y*Y + Z*Z))*q;

}

const Hypercomplex& Hypercomplex::log()

{

 double Length;

 Length = sqrt(X*X + Y*Y + Z*Z);

 Length = atan(Length/W);

 W = 0.0;

 X *= Length;

 Y *= Length;

 Z *= Length;

 return *this;

}

const Hypercomplex& Hypercomplex::u(Hypercomplex& q)

{

Hypercomplex u = q;

u.W = 0;

 return u;

}

const Hypercomplex Hypercomplex::pow(int k)

{

Hypercomplex b(1,0,0,0);

Hypercomplex a = *this;

 while (k)

{

 if (k & 1)

 {

  b *= a;

 }

 a *= a;

 k >>= 1;

}

 return b;

}

const Hypercomplex Hypercomplex::toHypercomplex(const double &d)

{

 double w, x, y, z;

w = d;

x = 0;

y = 0;

z = 0;

 return Hypercomplex(w,x,y,z);

}

const Hypercomplex Hypercomplex::toHypercomplex(const std::complex<double> &c)

{

 double w,x,y,z;

w = c.real();

x = c.imag();

y = 0;

z = 0;

 return Hypercomplex(w,x,y,z);

}

const Hypercomplex& Hypercomplex::Initialize(const double& w, const double& x, const double& y, const double& z)

{

W = w;

X = x;

Y = y;

Z = z;

 return *this;

}

void Hypercomplex::setW(const double &value)

{

W = value;

}

void Hypercomplex::setX(const double &value)

{

X = value;

}

void Hypercomplex::setY(const double &value)

{

Y = value;

}

void Hypercomplex::setZ(const double &value)

{

Z = value;

}

double Hypercomplex::getW()

{

 return W;

}

double Hypercomplex::getX()

{

 return X;

}

double Hypercomplex::getY()

{

 return Y;

}

double Hypercomplex::getZ()

{

 return Z;

}

Class DPF

#pragma once

#include <iostream>

#include <fstream>

#include <complex>

#include "Hypercomplex.h"

#include <conio.h>

#include <math.h>

using namespace std;

#ifndef PI

#define PI 3.14159265358979323846

#endif

class DPF

{

private:

 int N1, N2;

Hypercomplex w1, w2;

Hypercomplex** mas_w;

public:

DPF();

~DPF(void);

 double** loadMatrix(char name[]);

 void writeComplexMatrix(char name[], complex<double>** complex_matrix, int&);

 

 template<typename T>

T* binInverse(T* x, int& n)

{

 int dl = n/2;

 int st = n - 1;

 int j = 0;

 T s0 = 0;

 int k = 0;

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

 {

  if(i < j)

  {

   s0 = x[i];

   x[i] = x[j];

   x[j] = s0;

  }

  k = dl;

  while(k <= j)

  {

   j -= k;

   k >>= 1;

  }

  j += k;

 }

 return x;

}

 template<typename T>

T** binMtrInverse(T** matr, int& N)

{

 int size = 0;

  size = N;

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

  {

   binInverse(matr[i], size);

  }

  T* v = new T[size];

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

  {

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

   {

    v[j] = matr[j][i];

   }

   binInverse(v, size);

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

   {

    matr[j][i] = v[j];

   }

  }

  delete v;

 return matr;

}

Hypercomplex** matrInv(Hypercomplex** q, int& N);

 void genMatr(int N, char name[]);

Hypercomplex** _dpf(Hypercomplex** matrix);

Hypercomplex** Q44(Hypercomplex&, Hypercomplex&,  Hypercomplex&, Hypercomplex&, Hypercomplex&, Hypercomplex&,  

Hypercomplex&, Hypercomplex&,  Hypercomplex&, Hypercomplex&,  Hypercomplex&, Hypercomplex&,  Hypercomplex&, Hypercomplex&,  Hypercomplex&, Hypercomplex&);

Hypercomplex butterfly2(Hypercomplex&, Hypercomplex&, Hypercomplex&, Hypercomplex&, int&, int&, int&);

Hypercomplex Q_dft_basis4(Hypercomplex**, int&, const int&, const int&);

 int pos(int a);

Hypercomplex** HyperMatrix(const int&);

Hypercomplex** toHyperMatrix(double**, int&);

Hypercomplex** Q_ftfbs(Hypercomplex**, int &size);

Hypercomplex Q00(Hypercomplex**, int&, int&, int&);

Hypercomplex Qab(Hypercomplex**, int&, int&, int&, const int&, const int&);

 void shift(Hypercomplex**, int, int);

 void clcCoordinates(int&, int&, int, int);

 void clcMasW(int&);

 void clearMasW();

Hypercomplex** Q_algftf(Hypercomplex** , int&);

 

 void fft2(Hypercomplex**,int&);

 void fft4(Hypercomplex**,int&);

complex<double> toComplex(Hypercomplex&);

 const Hypercomplex automorphism(Hypercomplex&, const int&);

complex<double>** allocPartSpectr(Hypercomplex**, int&);

complex<double>** getComplexSpectr(Hypercomplex**, int&);

 void set(int& N);

 int get_N1();

 int get_N2();

 void setW1(int);

 void setW2(int);

};


Приложение Б. Принцип работы программы

На вход подаётся файл с числами, сгенерированный таким образом:

void DPF::genMatr(int N, char name[])

{

ofstream out(name);

out << N << ' ' << N << endl;

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

{

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

 {

  out << pow(3,double(k*j % 3)) << ' ';

 }

 out << endl;

}

 

out.close();

}.

Рассмотрим входной сигнал размером 4х4. Текстовый файл signal.txt имеет структуру:

4 4

1 1 1 1

1 3 9 1

1 9 3 1

1 1 1 1

Получаем результат в комплексной алгебре.

Результат БПФ по основанию 4, записанный в файл spectr_bpf.txt

(6,0) (10,2) (-2,0) (10,-2)

(10,2) (2,2) (-8,0) (4,0)

(-2,0) (-8,0) (2,0) (-8,0)

(10,-2) (4,0) (-8,0) (2,-2)

Результат ДПФ по основанию 4, записанный в файл spectr_dpf.txt

(36,0) (-10,10) (0,0) (-10,-10)

(-10,10) (0,-16) (6,6) (4,0)

(0,0) (6,6) (-12,0) (6,-6)

(-10,-10) (4,0) (6,-6) (0,16)


 

А также другие работы, которые могут Вас заинтересовать

1190. IP-адресация 78.5 KB
  Поиск и выбор IP-адреса. Центр управления сетями и общим доступом. Беспроводное сетевое соединение. Управление сетевыми подключениями.
1191. Фармакология периферической нервной системы 51 KB
  Структурно-функциональные особенности автономной (вегетативной) и соматической иннервации. Механизм передачи нервного импульса в синапсах автономной нервной системы. Эффекты возбуждения симпатических и парасимпатических нервов. Принципы фармакологической регуляции нейромедиаторных процессов. Классификация веществ, действующих на автономную нервную систему.
1192. Исследование ассортимента и качества конфет в розничной торговой сети «Рассвет» 346.34 KB
  Характеристика магазина «Рассвет» и основных показателей хозяйственной деятельности. Производители – источники поступления конфетных изделий в магазин «Рассвет». Приемка конфетных изделий по качеству и количеству в магазине «Рассвет». Анализ ассортимента конфетных изделий в магазине по видам. Результаты исследований ассортимента конфет в магазине по показателям: широта, полнота, устойчивость, новизна...
1193. Технико-экономические показатели автоматизации комплекса очистных сооружений 139.5 KB
  расчёт прогнозируемого объёма переработки загрязнённых вод на комплексе очистных сооружений. Расчет потребности капитальных вложений в основные фонды. Расчёт заработной платы. Расчёт затрат на потребление теплоэнергии (пара). Основные технико-экономические показатели автоматизации комплекса.
1194. Привод подвесного конвейера 1.04 MB
  Определение расчетной мощности электродвигателя. Выбор двигателя по каталогу. Определение передаточных чисел мощностей, частот вращения и крутящих моментов на валах привода. Расчет открытой прямозубой цилиндрической зубчатой передачи. Определение допускаемых напряжений для шестерни и колеса, при расчете на выносливость при изгибе.
1195. Организационно-экономическая часть создания системы автоматизированного проектирования 185.5 KB
  Технико-экономическое обоснование целесообразости проекта. Использование программно-аппаратных средств. Расчёт договорной цены разработки ячейки. Дополнительная заработная плата научного персонала. Календарный график работ по разработке блока.
1196. Газотурбинные установки. Машины и оборудование. Охрана труда. Строительные конструкции. 368.5 KB
  Устройство камер сгорания и теплообменных аппаратов ГТУ. Назначение, устройство и виды фильтров, используемых в гту. назначение и устройство глушителей, применяемых в компрессорах. Запорная арматура: назначение, устройство, принцип действия, примущества, недостатки, ремонтнопригодность. Устройство и принцип действия центробежного нагнетателя. Требования к проведению инструктажей по охране труда, их виды и сроки проведения.
1197. Приспособление для проведения механических испытаний 74.5 KB
  Схема сборки приспособления для проведения вибрационных испытаний. Универсальная оснастка плита для проведения вибрационных испытаний. Эта оснастка используется для проведения испытаний множества приборов.
1198. Производственный анализ СПК Октябрьский Волотовского района 63.5 KB
  Уставный капитал колхоза составляет 29531 тыс. руб. Учредителями хозяйства являются члены кооператива, которые объединили свои имущественные и земельные доли. Основным видом деятельности является сельское хозяйство, доля выручки от продажи сельскохозяйственной продукции на 2010 год составила 98 %.