Міністерство освіти і науки України
Національний університет „Львівська політехніка”
Кафедра ЕОМ
Звіт
з лабораторної роботи №2
на тему:
„Швидкі алгоритми обчислень дискретних тригонометричних перетворень.”
Підготував: ст. гр. СКС-12
Цигилик Л. О.
Прийняв: Кузьо М. М.
Львів 2007
Мета: Дослідити швидкі алгоритми дискретних тригонометричних перетворень і порівняти їх з алгоритмами безпосереднього обчислення тpигонометpних перетворень.
Теоретичні відомості:
Алгоритми швидкого перетворення Фур’є (ШПФ) можуть бути отримані за допомогою послідовного застосування операції розкладу одномірного масиву вхідних відліків сигналу на двохмірний. Ця операція здійснена тільки у випадку, коли N (довжина послідовності) є складним числом (N = N1 ( N2 ( ... ( Nj). Якщо N просте, його неможливо розкласти на прості співмножники; для такої послідовності алгоритмів ШПФ не існує. В більшості практичних випадків вхідну послідовність штучно продовжують додаванням нульових відліків до отримання N як складного числа.
Для характеристики розкладу використовують поняття “основа”, якщо всі співмножники однакові (N1 = N2 = ... = Nj) і “розщеплена основа”, якщо співмножники неоднакові.
Приклад. N = 64 = 2 ( 2 ( 2 ( 2 ( 2 ( 2 - основа 2.
N = 64 = 4 ( 4 ( 4 - основа 4.
N = 128 = 2 ( 4 ( 4 ( 4 - розщеплена основа 2-4.
Дискретне перетворення Фур’є
Дискретне перетворення Фур’є (ДПФ) кінцевої послідовності {x(n)}, 0 ( n ( N-1 визначається як
, k = 0, 1, ..., N-1 (1)
де .
Суть алгоритмів ШПФ в тому, що якщо N складне число і є степенем двійки (N=2m), то вихідна послідовність розбивається на дві коротші послідовності, ДПФ яких можуть бути скомбіновані таким чином, щоб утворилось ДПФ вихідної послідовності.
Методика побудови алгоритмів ШПФ наступна. Введемо дві (N/2) - точкові послідовності {x1(n)} і {x2(n)} з парних і непарних членів x(n) відповідно, x1(n) = x(2n) і x2(n) = x(2n + 1), n = 0, 1, ..., N/2-1. Тоді формулу розкладу можна записати так:
де X1(k) і X2(k) є N/2 - точкові ДПФ парних і непарних відліків вихідної послідовності.
Застосовуючи цю формулу розкладу X1(k) і X2(k) до пори, поки X1(k) і X2(k) не стануть двохточковим ДПФ, отримаємо алгоритм ШПФ з прорідженням в часі.
Для отримання іншої розповсюдженої форми алгоритмів ШПФ вихідну послідовність розбивають на дві (N/2) - точкові послідовності таким чином: перша складається з перших N/2 відліків, а друга з решти, x1(n) = x(n) і x2(n) = x(n + N/2), n = 0, 1, ..., N/2-1.
При такому розбитті N - точкове ДПФ можна записати так
де X1(k) і X2(k) є N/2 - точкові ДПФ першої та другої половини відліків вихідної послідовності.
Алгоритми ШПФ, що отримані з застосуванням цієї методики називаються алгоритмами з прорідженням по частоті.
Завдання до лабораторної роботи № 2.
Розробити процедуру дискретного тригонометричного перетворення вхідної послідовності x(n), розмірності N використовуючи алгоритм Algorithm. Формулу розкладу отримати за методом Кулі-Тьюкі. Algorithm : ШПФкі (швидке перетворення Фур’є за основою і) з прорідженням за часом (Т) чи частотою (F).
X(N) = ; N = 0, 1, ..., N-1. (1)
X(N) = ; N = 0, 1, ..., N-1. (2)
Варіант № 17
№ вар.
X(N)
Тип прорідження
Розмірність N
Основа
17
(1)
Т
64
2
Програмна реалізація завдання:
#include "stdafx.h"
#include "lab2csp.h"
#include <stdio.h>
#include <conio.h>
#include <malloc.h>
#include <stdlib.h>
#include <time.h>
// variant 17
#include <math.h>
#ifdef _DEBUG
#define new DEBUG_NEW
#undef THIS_FILE
static char THIS_FILE[] = __FILE__;
#endif
#define N 64
#define PI acos(-1)
// complex datatype
typedef struct{
long double rp;
long double ip;
} complex;
/////////////////////////////////////////////////////////////////////////////
// The one and only application object
//CWinApp theApp;
using namespace std;
//form input array
long double * FormF(long double *Xn)
{
//long double Xn[N];
for (int i=0;i<N;i++)
{
*(Xn+i) = 0;
for (int j=0;j<N/2;j++)
{
*(Xn+i) +=cos((2.0*PI*j*i/N));
};
*(Xn+i) *=1.0/N/2.0;
};
return Xn;
}
//DFT
complex *Dft(long double * Xn,complex *Xk)
{
// complex Xk[N];
for (int k=0;k<N;k++)
{
Xk[k].ip = 0;
Xk[k].rp = 0;
for (int n=0;n<N;n++)
{
Xk[k].rp += Xn[n]*cos((2.0*PI*k*n/(N*1.0)));
Xk[k].ip += Xn[n]*sin((2.0*PI*k*n/(N*1.0)));
};
Xk[k].ip*=-1;
};
return Xk;
}
// sort input array
long double * sortinp (long double * input, int k)
{
//for (int h=0; h<N; h++) printf("%e \n",input[h]);
long double *sorted= (long double *) malloc (sizeof (long double) * N);
int m=1;
for (int l=0;l<k;l++)
m*=2;
long double * arr1 = (long double *) malloc (sizeof (long double) * N / 2);
long double * arr2 = (long double *) malloc (sizeof (long double) * N / 2);
for (int i=0,j=0;i<N;j++)
{
arr1[j] = input[i];
//for (int l=0;l<k;l++)
i+=2;
};
for (int i=1,j=0;i<N;j++)
{
arr2[j] = input[i];
i+=2;
};
k++;
for (int i=0;i<N/2;i++)
{
sorted[i] = arr1[i];
};
for (int i=N/2;i<N;i++)
{
sorted[i] = arr2[i-N/2];
};
if (k<=5)
{
sorted=sortinp(sorted,k);
//arr2=sortinp(arr2,k);
}
free(arr1);
free(arr2);
return sorted;
}
// multiplying complex numbers
complex mulc(complex a,complex b)
{
complex c;
c.rp = a.rp*b.rp - a.ip*b.ip;
c.ip = a.ip*b.rp + b.ip*a.rp;
return c;
}
// adding complex numbers
complex addc(complex a,complex b)
{
complex c;
c.rp = a.rp+b.rp;
c.ip = a.ip+b.ip;
return c;
}
// subing complex numbers
complex subc(complex a,complex b)
{
complex c;
c.rp = a.rp-b.rp;
c.ip = a.ip-b.ip;
return c;
}
//FFT
complex * Fft(complex * Xn, int k)
{
complex * c =(complex*) calloc(64,4);//= new (512);
double n=exp((k-1)*log(2))/1;
k++;
complex w;
for (int i=0,j=0; i<N-n;)
{
w.ip = -1.0*(sin(2.0*PI*i/(N*1.0)));
w.rp = cos(2.0*PI*i/(N*1.0));
c[j] = addc(Xn[i], mulc(Xn[i+n],w));
c[j+1]= subc(Xn[i], mulc(Xn[i+n],w));
i+=n;
j+=2;
}
if (k<=6)
{
c=Fft(c,k);
}
return c;
}
int main(int argc, TCHAR* argv[], TCHAR* envp[])
{
int nRetCode = 0;
// initialize MFC and print and error on failure
if (!AfxWinInit(::GetModuleHandle(NULL), NULL, ::GetCommandLine(), 0))
// // TODO: change error code to suit your needs
{
cerr << _T("Fatal Error: MFC initialization failed") << endl;
nRetCode = 1;
}
else
{
long double *F=((long double *) malloc (sizeof (long double) * N)),*G=((long double *) malloc (sizeof (long double) * N));
FormF(F);
G=sortinp(F,1);
complex * fftc = ((complex *) malloc (sizeof (long double) * N*2)),* dftc = ((complex *) malloc (sizeof (long double) * N*2));
for (int j=0;j<N;j++)
{
fftc[j].ip = 0;
fftc[j].rp = G[j];
}
tm *tsd = (tm *) malloc (32),*tfd= (tm *) malloc (32),*tsf= (tm *) malloc (32),*tff= (tm *) malloc (32);
tm *td= (tm *) malloc (32),*tf= (tm *) malloc (32);
// tsd.GetCurrentTime();
_getsystime(tsd);
dftc=Dft(F,dftc);
_getsystime(tfd);
td->tm_sec = tsd->tm_sec - tfd->tm_sec;
// tsf.GetCurrentTime();
_getsystime(tff);
fftc = Fft(fftc,1);
_getsystime(tsf);
tf->tm_sec=tsf->tm_sec - tff->tm_sec;
FILE *fp = fopen("test","wb");
//printf("%f\n",acos(-1));
for (int i=0;i<N;i++)
{
fprintf(fp,"%02d %+08.8le %+08.8le*j %+08.8le %+08.8le*j\n",i+1,dftc[i].rp,dftc[i].ip,fftc[i].rp,fftc[i].ip);
};
fprintf(fp,"%d\t\t\t%d",td->tm_sec,tf->tm_sec);
fclose(fp);
getch();
};
return 0;
}
Висновок: у ході цієї лабораторної роботи, я дослідив швидкі алгоритми дискретних тригонометричних перетворень і порівняв їх з алгоритмами безпосереднього обчислення тpигонометpних перетворень.