1 // This file is part of the FXT library.
2 // Copyright (C) 2010, 2012 Joerg Arndt
3 // License: GNU General Public License version 3 or later,
4 // see the file COPYING.txt in the main directory.
5
6 #include "realfft/realfft.h"
7 #include "convolution/realconvolution.h"
8 #include "aux0/cmult.h"
9 #include "aux0/sumdiff.h"
10 #include "restrict.h"
11 #include "fxttypes.h" // ulong
12
13
14 void
fht_fft_convolution(double * restrict f,double * restrict g,ulong ldn)15 fht_fft_convolution(double * restrict f, double * restrict g, ulong ldn)
16 // (cyclic, real) convolution: g[] := f[] (*) g[]
17 // ldn := base-2 logarithm of the array length
18 {
19 fht_real_complex_fft(f, ldn);
20 fht_real_complex_fft(g, ldn);
21
22 fft_convolution_core1(f, g, ldn);
23
24 fht_complex_real_fft(g, ldn);
25 }
26 // -------------------------
27
28
29 void
split_radix_fft_convolution(double * restrict f,double * restrict g,ulong ldn)30 split_radix_fft_convolution(double * restrict f, double * restrict g, ulong ldn)
31 // (cyclic, real) convolution: g[] := f[] (*) g[]
32 // ldn := base-2 logarithm of the array length
33 // n = 2**ldn must be >=2
34 {
35 // if ( ldn==0 ) { g[0] *= f[0]; return; }
36
37 split_radix_real_complex_fft(f, ldn);
38 split_radix_real_complex_fft(g, ldn);
39
40 fft_convolution_core1(f, g, ldn);
41
42 split_radix_complex_real_fft(g, ldn);
43 }
44 // -------------------------
45
46
47 void
fht_fft_convolution0(double * restrict f,double * restrict g,ulong ldn)48 fht_fft_convolution0(double * restrict f, double * restrict g, ulong ldn)
49 // (linear, real) convolution: g[] := f[] (*) g[]
50 // ldn := base-2 logarithm of the array length
51 // input data must be zero padded:
52 // f[n/2] .. f[n-1] == 0 and g[n/2] .. g[n-1] == 0
53 // n = 2**ldn must be >=2
54 {
55 fht_real_complex_fft0(f, ldn);
56 fht_real_complex_fft0(g, ldn);
57
58 fft_convolution_core1(f, g, ldn);
59
60 fht_complex_real_fft(g, ldn);
61 }
62 // -------------------------
63
64
65 void
split_radix_fft_convolution0(double * restrict f,double * restrict g,ulong ldn)66 split_radix_fft_convolution0(double * restrict f, double * restrict g, ulong ldn)
67 // (linear, real) convolution: g[] := f[] (*) g[]
68 // ldn := base-2 logarithm of the array length
69 // input data must be zero padded:
70 // f[n/2] .. f[n-1] == 0 and g[n/2] .. g[n-1] == 0
71 // n = 2**ldn must be >=2
72 {
73 split_radix_real_complex_fft0(f, ldn);
74 split_radix_real_complex_fft0(g, ldn);
75
76 fft_convolution_core1(f, g, ldn);
77
78 split_radix_complex_real_fft(g, ldn);
79 }
80 // -------------------------
81
82
83
84 void
fft_convolution_core1(const double * restrict f,double * restrict g,ulong ldn,double v)85 fft_convolution_core1(const double * restrict f, double * restrict g,
86 ulong ldn, double v/*=0.0*/)
87 // Auxiliary routine for FFT based convolutions.
88 // Supply a value for v for a normalization factor != 1/n
89 {
90 const ulong n = (1UL<<ldn);
91
92 if ( v==0.0 ) v = 1.0/(double)n;
93
94 g[0] *= f[0] * v;
95 const ulong nh = n/2;
96 g[nh] *= f[nh] * v;
97
98 for (ulong i=1, j=n-1; i<j; ++i, --j)
99 {
100 cmult_n(f[i], f[j], g[i], g[j], v);
101 }
102 }
103 // -------------------------
104
105
106 void
fft_convolution_core2(const double * restrict f,double * restrict g,ulong ldn,double v)107 fft_convolution_core2(const double * restrict f, double * restrict g,
108 ulong ldn, double v/*=0.0*/)
109 // Auxiliary routine for FFT based convolutions.
110 // Supply a value for v for a normalization factor != 1/n
111 {
112 const ulong n = (1UL<<ldn);
113
114 if ( v==0.0 ) v = 1.0/(double)n;
115
116 g[0] *= f[0] * v;
117 const ulong nh = n/2;
118 g[nh] *= f[nh] * v;
119
120 for (ulong i=1, j=nh+1; i<nh; ++i, ++j)
121 {
122 cmult_n(f[i], f[j], g[i], g[j], v);
123 }
124 }
125 // -------------------------
126