1 /*
2 * Copyright (c) 2007 - 2015 Joseph Gaeddert
3 *
4 * Permission is hereby granted, free of charge, to any person obtaining a copy
5 * of this software and associated documentation files (the "Software"), to deal
6 * in the Software without restriction, including without limitation the rights
7 * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
8 * copies of the Software, and to permit persons to whom the Software is
9 * furnished to do so, subject to the following conditions:
10 *
11 * The above copyright notice and this permission notice shall be included in
12 * all copies or substantial portions of the Software.
13 *
14 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
15 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
16 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
17 * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
18 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
19 * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
20 * THE SOFTWARE.
21 */
22
23 //
24 // fft_r2r_1d.c : real-to-real methods (DCT/DST)
25 //
26
27 #include <stdio.h>
28 #include <stdlib.h>
29 #include <math.h>
30 #include "liquid.internal.h"
31
32 // create DCT/DST plan
33 // _nfft : FFT size
34 // _x : input array [size: _nfft x 1]
35 // _y : output array [size: _nfft x 1]
36 // _type : type (e.g. LIQUID_FFT_REDFT00)
37 // _method : fft method
FFT(_create_plan_r2r_1d)38 FFT(plan) FFT(_create_plan_r2r_1d)(unsigned int _nfft,
39 T * _x,
40 T * _y,
41 int _type,
42 int _flags)
43 {
44 // allocate plan and initialize all internal arrays to NULL
45 FFT(plan) q = (FFT(plan)) malloc(sizeof(struct FFT(plan_s)));
46
47 q->nfft = _nfft;
48 q->xr = _x;
49 q->yr = _y;
50 q->type = _type;
51 q->flags = _flags;
52
53 // TODO : use separate 'method' for real-to-real types
54 //q->method = LIQUID_FFT_METHOD_NONE;
55
56 switch (q->type) {
57 case LIQUID_FFT_REDFT00: q->execute = &FFT(_execute_REDFT00); break; // DCT-I
58 case LIQUID_FFT_REDFT10: q->execute = &FFT(_execute_REDFT10); break; // DCT-II
59 case LIQUID_FFT_REDFT01: q->execute = &FFT(_execute_REDFT01); break; // DCT-III
60 case LIQUID_FFT_REDFT11: q->execute = &FFT(_execute_REDFT11); break; // DCT-IV
61
62 case LIQUID_FFT_RODFT00: q->execute = &FFT(_execute_RODFT00); break; // DST-I
63 case LIQUID_FFT_RODFT10: q->execute = &FFT(_execute_RODFT10); break; // DST-II
64 case LIQUID_FFT_RODFT01: q->execute = &FFT(_execute_RODFT01); break; // DST-III
65 case LIQUID_FFT_RODFT11: q->execute = &FFT(_execute_RODFT11); break; // DST-IV
66 default:
67 fprintf(stderr,"error: fft_create_plan_r2r_1d(), invalid type, %d\n", q->type);
68 exit(1);
69 }
70
71 return q;
72 }
73
74 // destroy real-to-real transform plan
FFT(_destroy_plan_r2r_1d)75 void FFT(_destroy_plan_r2r_1d)(FFT(plan) _q)
76 {
77 // free main object memory
78 free(_q);
79 }
80
81 // print real-to-real transform plan
FFT(_print_plan_r2r_1d)82 void FFT(_print_plan_r2r_1d)(FFT(plan) _q)
83 {
84 printf("real-to-real transform...\n");
85 // TODO: print actual transform type
86 }
87
88 //
89 // DCT : Discrete Cosine Transforms
90 //
91
92 // DCT-I
FFT(_execute_REDFT00)93 void FFT(_execute_REDFT00)(FFT(plan) _q)
94 {
95 // ugly, slow method
96 unsigned int i,k;
97 float n_inv = 1.0f / (float)(_q->nfft-1);
98 float phi;
99 for (i=0; i<_q->nfft; i++) {
100 T x0 = _q->xr[0]; // first element
101 T xn = _q->xr[_q->nfft-1]; // last element
102 _q->yr[i] = 0.5f*( x0 + (i%2 ? -xn : xn));
103 for (k=1; k<_q->nfft-1; k++) {
104 phi = M_PI*n_inv*((float)k)*((float)i);
105 _q->yr[i] += _q->xr[k]*cosf(phi);
106 }
107
108 // compensate for discrepancy
109 _q->yr[i] *= 2.0f;
110 }
111 }
112
113 // DCT-II (regular 'dct')
FFT(_execute_REDFT10)114 void FFT(_execute_REDFT10)(FFT(plan) _q)
115 {
116 // ugly, slow method
117 unsigned int i,k;
118 float n_inv = 1.0f / (float)_q->nfft;
119 float phi;
120 for (i=0; i<_q->nfft; i++) {
121 _q->yr[i] = 0.0f;
122 for (k=0; k<_q->nfft; k++) {
123 phi = M_PI*n_inv*((float)k + 0.5f)*i;
124 _q->yr[i] += _q->xr[k]*cosf(phi);
125 }
126
127 // compensate for discrepancy
128 _q->yr[i] *= 2.0f;
129 }
130 }
131
132 // DCT-III (regular 'idct')
FFT(_execute_REDFT01)133 void FFT(_execute_REDFT01)(FFT(plan) _q)
134 {
135 // ugly, slow method
136 unsigned int i,k;
137 float n_inv = 1.0f / (float)_q->nfft;
138 float phi;
139 for (i=0; i<_q->nfft; i++) {
140 _q->yr[i] = _q->xr[0]*0.5f;
141 for (k=1; k<_q->nfft; k++) {
142 phi = M_PI*n_inv*((float)i + 0.5f)*k;
143 _q->yr[i] += _q->xr[k]*cosf(phi);
144 }
145
146 // compensate for discrepancy
147 _q->yr[i] *= 2.0f;
148 }
149 }
150
151 // DCT-IV
FFT(_execute_REDFT11)152 void FFT(_execute_REDFT11)(FFT(plan) _q)
153 {
154 // ugly, slow method
155 unsigned int i,k;
156 float n_inv = 1.0f / (float)(_q->nfft);
157 float phi;
158 for (i=0; i<_q->nfft; i++) {
159 _q->yr[i] = 0.0f;
160 for (k=0; k<_q->nfft; k++) {
161 phi = M_PI*n_inv*((float)k+0.5f)*((float)i+0.5f);
162 _q->yr[i] += _q->xr[k]*cosf(phi);
163 }
164
165 // compensate for discrepancy
166 _q->yr[i] *= 2.0f;
167 }
168 }
169
170 //
171 // DST : Discrete Sine Transforms
172 //
173
174 // DST-I
FFT(_execute_RODFT00)175 void FFT(_execute_RODFT00)(FFT(plan) _q)
176 {
177 // ugly, slow method
178 unsigned int i,k;
179 float n_inv = 1.0f / (float)(_q->nfft+1);
180 float phi;
181 for (i=0; i<_q->nfft; i++) {
182 _q->yr[i] = 0.0f;
183 for (k=0; k<_q->nfft; k++) {
184 phi = M_PI*n_inv*(float)((k+1)*(i+1));
185 _q->yr[i] += _q->xr[k]*sinf(phi);
186 }
187
188 // compensate for discrepancy
189 _q->yr[i] *= 2.0f;
190 }
191 }
192
193 // DST-II
FFT(_execute_RODFT10)194 void FFT(_execute_RODFT10)(FFT(plan) _q)
195 {
196 // ugly, slow method
197 unsigned int i,k;
198 float n_inv = 1.0f / (float)(_q->nfft);
199 float phi;
200 for (i=0; i<_q->nfft; i++) {
201 _q->yr[i] = 0.0f;
202 for (k=0; k<_q->nfft; k++) {
203 phi = M_PI*n_inv*((float)k+0.5f)*(i+1);
204 _q->yr[i] += _q->xr[k]*sinf(phi);
205 }
206
207 // compensate for discrepancy
208 _q->yr[i] *= 2.0f;
209 }
210 }
211
212 // DST-III
FFT(_execute_RODFT01)213 void FFT(_execute_RODFT01)(FFT(plan) _q)
214 {
215 // ugly, slow method
216 unsigned int i,k;
217 float n_inv = 1.0f / (float)(_q->nfft);
218 float phi;
219 for (i=0; i<_q->nfft; i++) {
220 _q->yr[i] = ((i%2)==0 ? 0.5f : -0.5f) * _q->xr[_q->nfft-1];
221 for (k=0; k<_q->nfft-1; k++) {
222 phi = M_PI*n_inv*((float)k+1)*((float)i+0.5f);
223 _q->yr[i] += _q->xr[k]*sinf(phi);
224 }
225
226 // compensate for discrepancy
227 _q->yr[i] *= 2.0f;
228 }
229 }
230
231 // DST-IV
FFT(_execute_RODFT11)232 void FFT(_execute_RODFT11)(FFT(plan) _q)
233 {
234 // ugly, slow method
235 unsigned int i,k;
236 float n_inv = 1.0f / (float)(_q->nfft);
237 float phi;
238 for (i=0; i<_q->nfft; i++) {
239 _q->yr[i] = 0.0f;
240 for (k=0; k<_q->nfft; k++) {
241 phi = M_PI*n_inv*((float)k+0.5f)*((float)i+0.5f);
242 _q->yr[i] += _q->xr[k]*sinf(phi);
243 }
244
245 // compensate for discrepancy
246 _q->yr[i] *= 2.0f;
247 }
248 }
249
250