libdspl-2.0
Библиотека алгоритмов цифровой обработки сигналов
fft.c
1 /*
2 * Copyright (c) 2015-2019 Sergey Bakhurin
3 * Digital Signal Processing Library [http://dsplib.org]
4 *
5 * This file is part of libdspl-2.0.
6 *
7 * is free software: you can redistribute it and/or modify
8 * it under the terms of the GNU Lesser General Public License as published by
9 * the Free Software Foundation, either version 3 of the License, or
10 * (at your option) any later version.
11 *
12 * DSPL is distributed in the hope that it will be useful,
13 * but WITHOUT ANY WARRANTY; without even the implied warranty of
14 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15 * GNU General Public License for more details.
16 *
17 * You should have received a copy of the GNU Lesser General Public License
18 * along with Foobar. If not, see <http://www.gnu.org/licenses/>.
19 */
20 
21 #include <stdlib.h>
22 #include <stdio.h>
23 #include <string.h>
24 #include "dspl.h"
25 #include "dspl_internal.h"
26 
27 /*******************************************************************************
28 COMPLEX vector IFFT
29 *******************************************************************************/
30 int DSPL_API ifft_cmplx(complex_t *x, int n, fft_t* pfft, complex_t* y)
31 {
32  int err, k;
33  double norm;
34 
35  if(!x || !pfft || !y)
36  return ERROR_PTR;
37  if(n<1)
38  return ERROR_SIZE;
39 
40 
41  err = fft_create(pfft, n);
42  if(err != RES_OK)
43  return err;
44 
45  memcpy(pfft->t1, x, n*sizeof(complex_t));
46  for(k = 0; k < n; k++)
47  IM(pfft->t1[k]) = -IM(pfft->t1[k]);
48 
49  err = fft_krn(pfft->t1, pfft->t0, pfft, n, 0);
50 
51  if(err!=RES_OK)
52  return err;
53 
54  norm = 1.0 / (double)n;
55  for(k = 0; k < n; k++)
56  {
57  RE(y[k]) = RE(pfft->t0[k])*norm;
58  IM(y[k]) = -IM(pfft->t0[k])*norm;
59  }
60  return RES_OK;
61 }
62 
63 
64 /*******************************************************************************
65 Real vector FFT
66 *******************************************************************************/
67 int DSPL_API fft(double *x, int n, fft_t* pfft, complex_t* y)
68 {
69  int err;
70 
71  if(!x || !pfft || !y)
72  return ERROR_PTR;
73  if(n<1)
74  return ERROR_SIZE;
75 
76 
77  err = fft_create(pfft, n);
78  if(err != RES_OK)
79  return err;
80 
81  re2cmplx(x, n, pfft->t1);
82 
83  return fft_krn(pfft->t1, y, pfft, n, 0);
84 }
85 
86 
87 
88 
89 /*******************************************************************************
90 COMPLEX vector FFT
91 *******************************************************************************/
92 int DSPL_API fft_cmplx(complex_t *x, int n, fft_t* pfft, complex_t* y)
93 {
94  int err;
95 
96  if(!x || !pfft || !y)
97  return ERROR_PTR;
98  if(n<1)
99  return ERROR_SIZE;
100 
101 
102  err = fft_create(pfft, n);
103  if(err != RES_OK)
104  return err;
105 
106  memcpy(pfft->t1, x, n*sizeof(complex_t));
107 
108  return fft_krn(pfft->t1, y, pfft, n, 0);
109 }
110 
111 
112 
113 /*******************************************************************************
114 composite FFT kernel
115 *******************************************************************************/
116 int fft_krn(complex_t* t0, complex_t* t1, fft_t* p, int n, int addr)
117 {
118  int n1, n2, k, m, i;
119  complex_t *pw = p->w+addr;
120  complex_t tmp;
121 
122  n1 = 1;
123  if(n%16== 0) { n1 = 16; goto label_size; }
124  if(n%7 == 0) { n1 = 7; goto label_size; }
125  if(n%5 == 0) { n1 = 5; goto label_size; }
126  if(n%4 == 0) { n1 = 4; goto label_size; }
127  if(n%3 == 0) { n1 = 3; goto label_size; }
128  if(n%2 == 0) { n1 = 2; goto label_size; }
129 
130 label_size:
131  if(n1 == 1)
132  {
133  for(k = 0; k < n; k++)
134  {
135  RE(t1[k]) = IM(t1[k]) = 0.0;
136  for(m = 0; m < n; m++)
137  {
138  i = (k*m) % n;
139  RE(tmp) = CMRE(t0[m], pw[i]);
140  IM(tmp) = CMIM(t0[m], pw[i]);
141  RE(t1[k]) += RE(tmp);
142  IM(t1[k]) += IM(tmp);
143  }
144  }
145  }
146  else
147  {
148  n2 = n / n1;
149 
150  if(n2>1)
151  {
152  memcpy(t1, t0, n*sizeof(complex_t));
153  transpose_cmplx(t1, n2, n1, t0);
154  }
155 
156  if(n1 == 16)
157  for(k = 0; k < n2; k++)
158  dft16(t0+16*k, t1+16*k);
159 
160  if(n1 == 7)
161  for(k = 0; k < n2; k++)
162  dft7(t0+7*k, t1+7*k);
163 
164  if(n1 == 5)
165  for(k = 0; k < n2; k++)
166  dft5(t0+5*k, t1+5*k);
167 
168  if(n1 == 4)
169  for(k = 0; k < n2; k++)
170  dft4(t0+4*k, t1+4*k);
171 
172  if(n1 == 3)
173  for(k = 0; k < n2; k++)
174  dft3(t0+3*k, t1+3*k);
175 
176  if(n1 == 2)
177  for(k = 0; k < n2; k++)
178  dft2(t0+2*k, t1+2*k);
179 
180  if(n2 > 1)
181  {
182 
183  for(k =0; k < n; k++)
184  {
185  RE(t0[k]) = CMRE(t1[k], pw[k]);
186  IM(t0[k]) = CMIM(t1[k], pw[k]);
187  }
188 
189  transpose_cmplx(t0, n1, n2, t1);
190 
191  for(k = 0; k < n1; k++)
192  {
193  fft_krn(t1+k*n2, t0+k*n2, p, n2, addr+n);
194  }
195  transpose_cmplx(t0, n2, n1, t1);
196  }
197  }
198  return RES_OK;
199 
200 }
201 
202 
203 
204 
205 /*******************************************************************************
206 FFT create for composite N
207 *******************************************************************************/
208 int DSPL_API fft_create(fft_t *pfft, int n)
209 {
210 
211  int n1, n2, addr, s, k, m, nw, err;
212  double phi;
213  s = n;
214  nw = addr = 0;
215 
216  if(pfft->n == n)
217  return RES_OK;
218 
219  while(s > 1)
220  {
221  n2 = 1;
222  if(s%16== 0) { n2 = 16; goto label_size; }
223  if(s%7 == 0) { n2 = 7; goto label_size; }
224  if(s%5 == 0) { n2 = 5; goto label_size; }
225  if(s%4 == 0) { n2 = 4; goto label_size; }
226  if(s%3 == 0) { n2 = 3; goto label_size; }
227  if(s%2 == 0) { n2 = 2; goto label_size; }
228 
229 
230 label_size:
231  if(n2 == 1)
232  {
233  if(s > FFT_COMPOSITE_MAX)
234  {
235  err = ERROR_FFT_SIZE;
236  goto error_proc;
237  }
238 
239  nw += s;
240  pfft->w = pfft->w ? (complex_t*) realloc(pfft->w, nw*sizeof(complex_t)):
241  (complex_t*) malloc( nw*sizeof(complex_t));
242  for(k = 0; k < s; k++)
243  {
244  phi = - M_2PI * (double)k / (double)s;
245  RE(pfft->w[addr]) = cos(phi);
246  IM(pfft->w[addr]) = sin(phi);
247  addr++;
248  }
249  s = 1;
250  }
251  else
252  {
253  n1 = s / n2;
254  nw += s;
255  pfft->w = pfft->w ? (complex_t*) realloc(pfft->w, nw*sizeof(complex_t)):
256  (complex_t*) malloc( nw*sizeof(complex_t));
257 
258  for(k = 0; k < n1; k++)
259  {
260  for(m = 0; m < n2; m++)
261  {
262  phi = - M_2PI * (double)(k*m) / (double)s;
263  RE(pfft->w[addr]) = cos(phi);
264  IM(pfft->w[addr]) = sin(phi);
265  addr++;
266  }
267  }
268  }
269  s /= n2;
270  }
271 
272  pfft->t0 = pfft->t0 ? (complex_t*) realloc(pfft->t0, n*sizeof(complex_t)):
273  (complex_t*) malloc( n*sizeof(complex_t));
274 
275  pfft->t1 = pfft->t1 ? (complex_t*) realloc(pfft->t1, n*sizeof(complex_t)):
276  (complex_t*) malloc( n*sizeof(complex_t));
277 
278  pfft->n = n;
279 
280  return RES_OK;
281 error_proc:
282  if(pfft->t0) free(pfft->t0);
283  if(pfft->t1) free(pfft->t1);
284  if(pfft->w) free(pfft->w);
285  pfft->n = 0;
286  return err;
287 }
288 
289 
290 
291 
292 
293 /*******************************************************************************
294 FFT free
295 *******************************************************************************/
296 void DSPL_API fft_free(fft_t *pfft)
297 {
298  if(!pfft)
299  return;
300  if(pfft->w)
301  free(pfft->w);
302  if(pfft->t0)
303  free(pfft->t0);
304  if(pfft->t1)
305  free(pfft->t1);
306  memset(pfft, 0, sizeof(fft_t));
307 }
308 
309 
310 
311 
312 /*******************************************************************************
313 FFT shifting
314 *******************************************************************************/
315 int DSPL_API fft_shift(double* x, int n, double* y)
316 {
317  int n2, r;
318  int k;
319  double tmp;
320  double *buf;
321 
322  if(!x || !y)
323  return ERROR_PTR;
324 
325  if(n<1)
326  return ERROR_SIZE;
327 
328  r = n%2;
329  if(!r)
330  {
331  n2 = n>>1;
332  for(k = 0; k < n2; k++)
333  {
334  tmp = x[k];
335  y[k] = x[k+n2];
336  y[k+n2] = tmp;
337  }
338  }
339  else
340  {
341  n2 = (n-1) >> 1;
342  buf = (double*) malloc(n2*sizeof(double));
343  memcpy(buf, x, n2*sizeof(double));
344  memcpy(y, x+n2, (n2+1)*sizeof(double));
345  memcpy(y+n2+1, buf, n2*sizeof(double));
346  free(buf);
347  }
348  return RES_OK;
349 }
350 
351 
352 
353 /*******************************************************************************
354 FFT shifting for complex vector
355 *******************************************************************************/
356 int DSPL_API fft_shift_cmplx(complex_t* x, int n, complex_t* y)
357 {
358  int n2, r;
359  int k;
360  complex_t tmp;
361  complex_t *buf;
362 
363  if(!x || !y)
364  return ERROR_PTR;
365 
366  if(n<1)
367  return ERROR_SIZE;
368 
369  r = n%2;
370  if(!r)
371  {
372  n2 = n>>1;
373  for(k = 0; k < n2; k++)
374  {
375  RE(tmp) = RE(x[k]);
376  IM(tmp) = IM(x[k]);
377 
378  RE(y[k]) = RE(x[k+n2]);
379  IM(y[k]) = IM(x[k+n2]);
380 
381  RE(y[k+n2]) = RE(tmp);
382  IM(y[k+n2]) = IM(tmp);
383  }
384  }
385  else
386  {
387  n2 = (n-1) >> 1;
388  buf = (complex_t*) malloc(n2*sizeof(complex_t));
389  memcpy(buf, x, n2*sizeof(complex_t));
390  memcpy(y, x+n2, (n2+1)*sizeof(complex_t));
391  memcpy(y+n2+1, buf, n2*sizeof(complex_t));
392  free(buf);
393  }
394  return RES_OK;
395 }
396 
complex_t * w
Definition: dspl.h:47
#define ERROR_SIZE
Ошибка при передаче размера массива.
Definition: dspl.h:125
Структура данных объекта быстрого преобразования Фурье
Definition: dspl.h:45
double complex_t[2]
Описание комплексного типа данных.
Definition: dspl.h:41
#define ERROR_PTR
Ошибка указателя.
Definition: dspl.h:118
complex_t * t1
Definition: dspl.h:49
#define IM(x)
Макрос определяющий мнимую часть комплексного числа.
Definition: dspl.h:65
#define RE(x)
Макрос определяющий реальную часть комплексного числа.
Definition: dspl.h:64
int fft_shift(double *x, int n, double *y)
Перестановка спектральных отсчетов дискретного преобразования Фурье
Definition: fft.c:315
complex_t * t0
Definition: dspl.h:48
int re2cmplx(double *x, int n, complex_t *y)
Преобразование массива вещественных данных в массив комплексных данных.
Definition: complex.c:82
int n
Definition: dspl.h:50
void fft_free(fft_t *pfft)
Очистить структуру fft_t алгоритма БПФ
Definition: fft.c:296
#define ERROR_FFT_SIZE
Неверно задан размер БПФ.
Definition: dspl.h:94
#define RES_OK
Функция завершилась корректно. Ошибки отсутствуют.
Definition: dspl.h:81
int fft_create(fft_t *pfft, int n)
Заполнение структуры fft_t для алгоритма БПФ
Definition: fft.c:208