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
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
|
/**********************************************************************
SLOW FOURIER TRANSFORM
wn_sft_vect(vector,len_i)
wn_inverse_sft_vect(vector,len_i)
**********************************************************************/
/****************************************************************************
COPYRIGHT NOTICE:
The source code in this file is provided free of charge
to the author's consulting clients. It is in the
public domain and therefore may be used by anybody for
any purpose.
AUTHOR:
Will Naylor
****************************************************************************/
#include <math.h>
#include "wnlib.h"
#include "wnasrt.h"
#include "wnmem.h"
#include "wncplx.h"
#include "wnfft.h"
local wn_cplx *roots_of_unity,*in_copy;
void wn_inverse_sft_vect(wn_cplx vector[],int len_i)
{
int i,j,exponent;
struct wn_cplx_struct prod_struct;
wn_cplx sum,prod;
double norm_factor;
wn_gpmake("no_free");
wn_cplx_make_vect(&roots_of_unity,len_i);
wn_cplx_make_vect(&in_copy,len_i);
/* compute roots of unity */
for(i=0;i<len_i;++i)
{
wn_polar_to_cplx(roots_of_unity[i],
1.0,2.0*M_PI*((double)i)/((double)len_i));
}
wn_cplx_copy_vect(in_copy,vector,len_i);
norm_factor = 1.0/sqrt((double)len_i);
prod = &prod_struct;
for(i=0;i<len_i;i++)
{
sum = vector[i];
sum->real = 0.0;
sum->imag = 0.0;
exponent = 0;
for(j=0;j<len_i;++j)
{ /* this is the inner loop */
wn_cplx_multiply(prod,in_copy[j],roots_of_unity[exponent]);
sum->real += prod->real;
sum->imag += prod->imag;
exponent +=i;
if(exponent >= len_i)
{
exponent -= len_i;
}
}
sum->real *= norm_factor;
sum->imag *= norm_factor;
}
wn_gpfree();
}
void wn_sft_vect(wn_cplx vector[],int len_i)
{
wn_cplx_conjugate_vect(vector,len_i);
wn_inverse_sft_vect(vector,len_i);
wn_cplx_conjugate_vect(vector,len_i);
}
|
