-
Notifications
You must be signed in to change notification settings - Fork 0
/
dft.cpp
68 lines (64 loc) · 2.14 KB
/
dft.cpp
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
/*
Äèñêðåòíîå ïðåîáðàçîâàíèå ôóðüå
*/
//---------------------------------------------------------------------------
/*âñå çàèíêëóäåííûå õèäåðû äî #pragma hdrstop áóäóò ïðåêîìïèëèðîâàíû è çàêåøåíû,
åñëè ãäå åùå â C/CPP ôàéëàõ êîìïèëÿòîð âñòðåòèò òó æå ïîñëåäîâàòåëüíîñòü
õèäåðîâ äî #pragma hdrstop - òî çàþçàåòñÿ óæå ñêîìïèëåíûé êåø.*/
#pragma hdrstop
#include "dft.h"
#pragma package(smart_init)
//---------------------------------------------------------------------------
DSP::DSP(double &T, unsigned int N){
this->T = T;
this->N = N;
this->w = 2*M_PI/(N*T);
this->q = w*T;
this->sample = new double[this->N];
this->M = 0;
this->oldest = 0;
}
//---------------------------------------------------------------------------
DSP::~DSP(){
delete [] this->sample;
}
//---------------------------------------------------------------------------
Vector FFT::calc(double &new_samlpe){
double x=0, y=0;
double q_m; //ñåðåäèíà âûáîðêè â ðàäèàíàõ
if(M < N){
sample[M] = new_samlpe;
M++; //÷èñëî ýëåìåíòîâ â î÷åðåäè
oldest = 0;
q_m = (M-1)*q/2;
}
else{
sample[oldest] = new_samlpe;
q_m = 0;
oldest = (++oldest)%N;
}
for(unsigned int n=0, i=oldest; n < M; n++, i=(++i)%N)
x += sample[i]*sin(q*n-q_m);
x = 2*x/N;
//êîýôèöèåíò äëÿ ôèëüòðàöèè ïîñòîÿííîé ñîñòàâëÿþùåé
double A_d = sin(q*M/2)/(M*sin(q/2));
//ïîïðàâî÷íûé êîýôôèöèåíò èç-çà ïîäàâëåíèÿ ïîñòîÿííîé ñîñòàâëÿþùåé
double tmp = 1.0*M/N + sin(q*M)/sin(q)/N;
double C_gd = 1.0/(1.0 - 2.0/N/M*pow(sin(q*M/2.0),2)/pow(sin(q/2.0),2)/tmp );
for(unsigned int n=0, i=oldest; n < M; n++, i=(++i)%N)
y += sample[i]*(cos(q*n-q_m)-A_d); //ñ ó÷åòîì ïîäàâëåíèÿ ïîñòîÿííîé ñîñòîâëÿþùåé;
y = 2*y/N;
y = y * C_gd;
// Ïîïðàâî÷íûå êîýôèöèåíòû, ò.ê. M ìîæåò áûòü ìåíüøå N
if( M<N && M>1 ){
double C_h = N/(M - sin(2*M_PI*M/N)/sin(q));
double C_g = N/(M + sin(2*M_PI*M/N)/sin(q));
x = x * C_h;
y = y * C_g;
}
Vector ret(x,y);
// ñ ó÷åòîì ïðèâåäåíèÿ óãëà
if(M<=N && x!=0) // <= ò.ê. ýòî ââåðõó áûëî M++, ò.å. ýòî âñå åùå íåïîëíîïåðèîäíûé
{double a=-q_m+M*q; ret.rotate(a);}
return ret;
}