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

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

{

Real[i*N+j] /= temp;

Image[i*N+j] /= temp;

}

delete[] Cos;

delete[] Sin;

}

//БПФ по строчкам и столбцам, строчки и столбцы уже переставлены

//для главного процесса передаются данные, для остальных Real = NULL, Image = NULL

void Fure2InLine(int Height, int log2Height, int Width, int log2Width, float *RealLine, float *ImageLine, bool flag, int rank, int size)

{

int i, j, oper, HeightPerProcess, WidthPerProcess;

int PositionReal, PositionImage, SizeOfColumn, SizeOfColumnAsLine, tag;

float *Real_Array_Height, *Image_Array_Height, *Real_Array_Width, *Image_Array_Width,

*Real_Array_Height_Packed, *Image_Array_Height_Packed, *ptrReal = NULL, *ptrImage = NULL;

if (flag)

oper = 1;

else

oper = -1;

HeightPerProcess = Height/size; //целое число - проверяется в главной функции

WidthPerProcess = Width/size;

Real_Array_Width = new float[Width];

Image_Array_Width = new float[Width];

float *Cos_Height, *Sin_Height, *Cos_Width, *Sin_Width, arg, TempSin, TempCos;//, temp;

Cos_Height = new float[(Height+1)/2]; //проверка при длине массива в 1

Sin_Height = new float [(Height+1)/2]; // [0,Pi]

arg = 2*M_PI/(float)Height;

TempSin = sin(arg);

TempCos = cos(arg);

Cos_Height[0] = 1; //cos(0)=1

Sin_Height[0] = 0; //sin(0) = 0

if(Height != 2)

{

Cos_Height[1] = TempCos;

Sin_Height[1] = TempSin;

}

for(i=2;i<Height/2;i++)

{

Cos_Height[i] = Cos_Height[i-1]*TempCos - Sin_Height[i-1]*TempSin;

Sin_Height[i] = Sin_Height[i-1]*TempCos + Cos_Height[i-1]*TempSin;

}

Cos_Width = new float [(Width+1)/2]; //проверка при длине массива в 1

Sin_Width = new float [(Width+1)/2]; // [0,Pi]

arg = 2*M_PI/(float)Width;

TempSin = sin(arg);

TempCos = cos(arg);

Cos_Width[0] = 1; //cos(0)=1

Sin_Width[0] = 0; //sin(0) = 0

if(Width!=2)

{

Cos_Width[1] = TempCos;

Sin_Width[1] = TempSin;

}

for(i=2;i<Width/2;i++) //считаем текущую комплексную экспоненту поворотом предыдущей на угол

{

Cos_Width[i] = Cos_Width[i-1]*TempCos - Sin_Width[i-1]*TempSin;

Sin_Width[i] = Sin_Width[i-1]*TempCos + Cos_Width[i-1]*TempSin;

}

ptrReal = RealLine; //ставим указатель для пересылки на начало данных

ptrImage = ImageLine;

for(i=0;i<HeightPerProcess;i++) //по строчкам

{

MPI_Barrier(MPI_COMM_WORLD);

//отправляем по строчке данных

MPI_Scatter(ptrReal, Width,MPI_FLOAT,Real_Array_Width, Width,MPI_FLOAT,0,MPI_COMM_WORLD); MPI_Scatter(ptrImage,Width,MPI_FLOAT,Image_Array_Width,Width,MPI_FLOAT,0,MPI_COMM_WORLD);

//считам БПФ над каждой строчкой Fure_line_with_trig(Width,log2Width,Real_Array_Width,Image_Array_Width,Cos_Width,Sin_Width,oper); //считам БПФ

MPI_Barrier(MPI_COMM_WORLD);

//принимаем по строчке данных

MPI_Gather(Real_Array_Width, Width,MPI_FLOAT,ptrReal, Width,MPI_FLOAT,0,MPI_COMM_WORLD); MPI_Gather(Image_Array_Width,Width,MPI_FLOAT,ptrImage,Width,MPI_FLOAT,0,MPI_COMM_WORLD);

if(rank == 0)

{

ptrReal+=Width; //сдвигаем указатель на следующую строку

ptrImage+=Width;

}

}

//транспонируем данные

int temp;

if(rank == 0)

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

for(j=i+1;j<Width;j++)

{

temp = RealLine[i*Height+j];

RealLine[i*Height+j] = RealLine[j*Width+i];

RealLine[j*Width+i] = temp;

temp = ImageLine[i*Height+j];

ImageLine[i*Height+j] = ImageLine[j*Width+i];

ImageLine[j*Width+i] = temp;

}

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

ptrReal = RealLine; //ставим указатель для пересылки на начало данных

ptrImage = ImageLine;

for(i=0;i<HeightPerProcess;i++) //по строчкам (столбцам)

{

MPI_Barrier(MPI_COMM_WORLD);

//отправляем по строчке данных

MPI_Scatter(ptrReal, Width,MPI_FLOAT,Real_Array_Width, Width,MPI_FLOAT,0,MPI_COMM_WORLD); MPI_Scatter(ptrImage,Width,MPI_FLOAT,Image_Array_Width,Width,MPI_FLOAT,0,MPI_COMM_WORLD);

//считам БПФ над каждой строчкой

Fure_line_with_trig(Width,log2Width,Real_Array_Width,Image_Array_Width,Cos_Width,Sin_Width,oper); //считам БПФ

MPI_Barrier(MPI_COMM_WORLD);

//принимаем по строчке данных

MPI_Gather(Real_Array_Width, Width,MPI_FLOAT,ptrReal, Width,MPI_FLOAT,0,MPI_COMM_WORLD);

MPI_Gather(Image_Array_Width,Width,MPI_FLOAT,ptrImage,Width,MPI_FLOAT,0,MPI_COMM_WORLD);

if(rank == 0)

{

ptrReal+=Width; //сдвигаем указатель на следующую строку

ptrImage+=Width;

}

}

delete[] Cos_Height;

delete[] Sin_Height;

delete[] Cos_Width;

delete[] Sin_Width;

delete[] Real_Array_Width;

delete[] Image_Array_Width;

}

#endif //FOURIER_CPP

Файл main.cpp - ввод и сохранение данных, замер времени, вызов преобразования Фурье.

#include <stdio.h>

#include <conio.h>

#include <math.h>

#include <time.h>

#include <iostream>

#include <mpi.h>

using namespace std;

const double Pi = 3.1415926535897932384626433832795;

//БПФ по строчкам и столбцам

void Fure2(int Height, int log2Height, int Width, int log2Width, float **Real, float **Image, bool flag, int rank, int size);

//БПФ по строчкам и столбцам

void Fure2InLine(int Height, int log2Height, int Width, int log2Width, float *RealLine, float *ImageLine, bool flag, int rank, int size);

//модифицированный алгоритм Кули-Тьюки

void Fourier_Kyli_Tuki (int N, int s, float *Real, float *Image, bool flag, int rank, int size);

inline int IsLog2(int N) //return -1 if N is not degree of 2

{

int i=1,s=0;

for(i=1;i<N;i*=2)

s++;

if(i == N)

return s;

else

return -1;

}

int revs(int i, int s)

{

int k;

if(i==0 && s==0)

return 0;

else

{

if(int(i/2)==i/2.0)

k=revs(int(i/2), s-1);

else

k=revs(int(i/2), s-1)+pow(2.0, s-1);

return k;

}

}

bool Save(char *FileName, int Height, int Width, float **Real, float **Image)

{

FILE *out;

out = fopen(FileName,"w");

if(out == NULL)

{

cout << "Error! File " << FileName << " cann't be open for write!" << endl;

return false;

}

int i,j;

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

{

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

if(Image[i][j] == 0)

fprintf(out,"%.2f ",Real[i][j]);

else

if(Image[i][j] < 0)

fprintf(out,"%.2f%.2fi ",Real[i][j],Image[i][j]);

else

fprintf(out,"%.2f+%.2fi ",Real[i][j],Image[i][j]);

fprintf(out,"\n");

}

fclose(out);

return true;

}

bool Save(char *FileName, int Height, int Width, float *Real, float *Image)

{

FILE *out;

out = fopen(FileName,"w");

if(out == NULL)

{

cout << "Error! File " << FileName << " cann't be open for write!" << endl;

return false;

}

int i,j,k;

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

{

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

if(Image[k] == 0)

fprintf(out,"%.2f ",Real[k]);

else

if(Image[k] < 0)

fprintf(out,"%.2f%.2fi ",Real[k],Image[k]);

else

fprintf(out,"%.2f+%.2fi ",Real[k],Image[k]);

fprintf(out,"\n");

}

fclose(out);

return true;

}

int main(int argc, char* argv[])

{

int size,rank;

MPI_Init(&argc,&argv);

MPI_Comm_size(MPI_COMM_WORLD, &size); // number of processes involved in run

MPI_Comm_rank(MPI_COMM_WORLD, &rank); // my process id: 0 <= rank < size

if(IsLog2(size)<0)

{

cout << "Error! " << size << " is not degree of 2!" << endl;

MPI_Finalize();

return -1;

}

if(rank == 0)

cout << "Number of process is " << size << endl;

FILE *in,*out,*f;

int Height = 0, Width = 0,

tempHeight, tempWidth,

Log2Height = 0, Log2Width = 0,

i, j;

_int64 k;

char FileNameInput[20] = "input.txt",

FileNameParam[20] = "param.txt",

FileNameOutput[20] = "output.txt",

FileNameKT[20] = "KT.txt",

FileNameLC[20] = "LC.txt";

float **Real = NULL, **Image = NULL, **RealData = NULL,

*RealLine = NULL, *ImageLine = NULL;

int Choise;

clock_t begin,end,WorkTime,CreateMemoryTime,LoadDataTime,SaveDataTime,GenerationDataTime,

ReversionTime, FFT2KTTime, FFT2LCTime,SaveFFT2KT,SaveFFT2LC;

float second = 1.0 ; //for calc second, 1000 for Win, 1000000 for Linux

int RANGE_MIN = 0, RANGE_MAX = 255;

int *IndexHeight=NULL,*IndexWidth=NULL;

if(rank == 0)

{

cout << "1 - Create data2" << endl;

cout << "2 - Load data2" << endl;

cin >> Choise;

switch(Choise)

{

case 1:

//Create data2

cout << "Input Height: ";

cin >> Height;

Log2Height = IsLog2(Height);

if(Log2Height < 0)

{

cout << "Error in height!" << endl;

return -1;

}

cout << "Input Width: ";

cin >> Width;

Log2Width= IsLog2(Width);

if(Log2Width < 0)

{

cout << "Error in Width!" << endl;

return -1;

}

f = fopen(FileNameParam,"w");

//f = fopen("param.txt","w");

if(f == NULL)

{

cout << "Error in open file " << FileNameParam << " for write parametrs" << endl;

return -1;

}

fprintf(f,"%d %d\n",Height,Width);

fclose(f);

begin = clock();

RealData = new float*[Height];

Real = new float*[Height];

Image = new float*[Height];

//FourierReal = new float*[Height];

//FourierImage = new float*[Height];

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

{

RealData[i] = new float[Width];

Real[i] = new float[Width];

Image[i] = new float[Width];

//FourierReal[i] = new float[Width];

//FourierImage[i] = new float[Width];

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

{

Real[i][j] = 0;

Image[i][j] = 0;

//FourierReal[i][j] = 0;

//FourierImage[i][j] = 0;

}

}

RealLine = new float[Height*Width];

ImageLine = new float[Height*Width];

for(k=0;k<Height*Width;k++)

{

RealLine[k] = 0;

ImageLine[k] = 0;

}

end = clock();

CreateMemoryTime = (end - begin)/second;

cout << "CreateMemoryTime = " << CreateMemoryTime << endl; begin = clock();

srand( (unsigned)time( NULL ) );

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

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

RealData[i][j] = (((float) rand() / (float) RAND_MAX) * RANGE_MAX + RANGE_MIN);

end = clock();

GenerationDataTime = (end - begin)/second;

cout << "GenerationDataTime = " << GenerationDataTime << endl;

in = fopen(FileNameInput,"w");

if(in == NULL)

{

cout << "Error! File " << FileNameInput << " cann't be open for write data!" << endl;

return -1;

}

begin = clock();

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

{

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

fprintf(in,"%d ",(int)RealData[i][j]);

fprintf(in,"\n");

}

end = clock();

fclose(in);

SaveDataTime = (end - begin)/second;

cout << "SaveDataTime = " << SaveDataTime << endl;

break;

case 2:

//Load data2

f = fopen(FileNameParam,"r");

if(f == NULL)

{

cout << "Error! File " << FileNameParam<< " cann't be open for read data!" << endl;

//need delete memory in ideal case

return -1;

}

fscanf(f,"%d%d",&Height,&Width);

fclose(f);

Log2Height = IsLog2(Height);

Log2Width = IsLog2(Width);

if(Log2Height<0 && Log2Width<0)

{

cout << "Error! Incorrect parametrs!" << endl;

MPI_Finalize();

return -1;

}

begin = clock();

RealData = new float*[Height];

Real = new float*[Height];

Image = new float*[Height];

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

{

RealData[i] = new float[Width];

Real[i] = new float[Width];

Image[i] = new float[Width];

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

{

Real[i][j] = 0;

Image[i][j] = 0;

}

}

RealLine = new float[Height*Width];

ImageLine = new float[Height*Width];

for(k=0;k<Height*Width;k++)

{

RealLine[k] = 0;

ImageLine[k] = 0;

}

end = clock();

CreateMemoryTime = (end - begin)/second;

cout << "CreateMemoryTime = " << CreateMemoryTime << endl;

in = fopen(FileNameInput,"r");

if(in == NULL)

{

cout << "Error! File " << FileNameInput << " cann't be open for read data!" << endl;

//need delete memory in ideal case

return -1;

}

begin = clock();

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

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

fscanf(in,"%f",&RealData[i][j]);

end = clock();

LoadDataTime = (end - begin)/second;

cout << "LoadDataTime = " << LoadDataTime << endl;

fclose(in);

break;

default:

cout <<"Error in Choise!\n" << endl;

MPI_Finalize();

return -1;

}

cout << "Height = " << Height << ", Width = " << Width << endl;

IndexHeight = new int[Height];

IndexWidth = new int[Width];

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

IndexHeight[i] = revs(i,Log2Height);

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

IndexWidth[j] = revs(j,Log2Width);

begin = clock();

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

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

{

Real[i][j] = RealData[IndexHeight[i]][IndexWidth[j]];

RealLine[k] = RealData[IndexHeight[i]][IndexWidth[j]];

}

end = clock();

ReversionTime = (end - begin)/second;

} //rank == 0

MPI_Barrier(MPI_COMM_WORLD);

MPI_Bcast(&Height,1,MPI_INT,0,MPI_COMM_WORLD);

MPI_Barrier(MPI_COMM_WORLD);

MPI_Bcast(&Width,1,MPI_INT,0,MPI_COMM_WORLD);

MPI_Barrier(MPI_COMM_WORLD);

MPI_Bcast(&Log2Height,1,MPI_INT,0,MPI_COMM_WORLD);

MPI_Barrier(MPI_COMM_WORLD);

MPI_Bcast(&Log2Width,1,MPI_INT,0,MPI_COMM_WORLD);

if(rank == 0)

{

cout << "FFT2 Line and Columns" << endl;

begin = clock();

}

Fure2InLine(Height,Log2Height,Width,Log2Width,RealLine,ImageLine,true,rank,size);

if(rank == 0)

{

end = clock();

FFT2LCTime = (end - begin)/second;

cout << "FFT2LCTime = " << FFT2LCTime << endl;

}

if(rank == 0)

{

begin = clock();

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

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

{

RealLine[k] = RealData[IndexHeight[i]][IndexWidth[j]];

ImageLine[k] = 0;

}

end = clock();

ReversionTime = (end - begin)/second;

}

if(rank == 0)

{

cout << "FFT2 Kyli-Tuki" << endl;

begin = clock();

}

Fourier_Kyli_Tuki(Height,Log2Height,RealLine,ImageLine,true, rank, size);

if(rank == 0)

{

end = clock();

FFT2KTTime = (end - begin)/second;

cout << "FFT2KTTime = " << FFT2KTTime << endl;

}

MPI_Barrier(MPI_COMM_WORLD);

cout << rank << ": End of program" << endl;

MPI_Finalize();

return 1;

}

Файл param.txt - хранит два числа: размеры сигнала, на данный момент программа работает только с «квадратным» сигналом, сторона которого - степень 2.

Файл input.txt - содержит сам сигнал.

42 61 207 251

115 164 26 174

107 77 73 33

156 124 86 150

Файл KT.txt - содержит результат двумерного преобразования Фурье по аналогу Кули-Тьюки

588.51 -83.56-51.44i 24.99 -83.56+51.44i

-92.65+9.86i -73.75+46.17i 11.14+56.38i 66.46+2.84i

-100.01 8.19-59.84i -27.22 8.19+59.84i

-92.65-9.86i 66.46-2.84i 11.14-56.38i -73.75-46.17i



Файл LC.txt - содержит результат двумерного преобразования Фурье по строчкам и столбцам

588.51 -92.65+9.86i -100.01 -92.65-9.86i

-83.40-50.96i -73.59+46.66i 8.36-59.35i 66.63-2.35i

25.10 11.00+56.38i -27.10 11.00-56.38i

-83.19+50.95i 66.45+2.69i 8.19+59.05i -73.45-46.69i

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