libdspl-2.0
Digital Signal Processing Algorithm Library
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 <float.h>
25 #include "dspl.h"
26 #include "dspl_internal.h"
27 
28 
29 #ifdef DOXYGEN_ENGLISH
30 
103 #endif
104 #ifdef DOXYGEN_RUSSIAN
105 
178 #endif
179 int DSPL_API ifft_cmplx(complex_t *x, int n, fft_t* pfft, complex_t* y)
180 {
181  int err, k;
182  double norm;
183 
184  if(!x || !pfft || !y)
185  return ERROR_PTR;
186  if(n<1)
187  return ERROR_SIZE;
188 
189 
190  err = fft_create(pfft, n);
191  if(err != RES_OK)
192  return err;
193 
194  memcpy(pfft->t1, x, n*sizeof(complex_t));
195  for(k = 0; k < n; k++)
196  IM(pfft->t1[k]) = -IM(pfft->t1[k]);
197 
198  err = fft_krn(pfft->t1, pfft->t0, pfft, n, 0);
199 
200  if(err!=RES_OK)
201  return err;
202 
203  norm = 1.0 / (double)n;
204  for(k = 0; k < n; k++)
205  {
206  RE(y[k]) = RE(pfft->t0[k])*norm;
207  IM(y[k]) = -IM(pfft->t0[k])*norm;
208  }
209  return RES_OK;
210 }
211 
212 
213 #ifdef DOXYGEN_ENGLISH
214 
283 #endif
284 #ifdef DOXYGEN_RUSSIAN
285 
355 #endif
356 int DSPL_API fft(double* x, int n, fft_t* pfft, complex_t* y)
357 {
358  int err;
359 
360  if(!x || !pfft || !y)
361  return ERROR_PTR;
362  if(n<1)
363  return ERROR_SIZE;
364 
365 
366  err = fft_create(pfft, n);
367  if(err != RES_OK)
368  return err;
369 
370  re2cmplx(x, n, pfft->t1);
371 
372  return fft_krn(pfft->t1, y, pfft, n, 0);
373 }
374 
375 
376 
377 #ifdef DOXYGEN_ENGLISH
378 
379 #endif
380 #ifdef DOXYGEN_RUSSIAN
381 
382 #endif
383 int DSPL_API fft_abs(double* x, int n, fft_t* pfft,
384  double fs, int flag,
385  double* mag, double* freq)
386 {
387  int k, err = RES_OK;
388  complex_t *X = NULL;
389  if(!x || !pfft)
390  return ERROR_PTR;
391 
392  if(n<1)
393  return ERROR_SIZE;
394 
395  if(mag)
396  {
397  X = (complex_t*)malloc(n*sizeof(complex_t));
398  err = fft(x, n, pfft, X);
399  if(err!=RES_OK)
400  goto error_proc;
401 
402 
403  for(k = 0; k < n; k++)
404  mag[k] = ABS(X[k]);
405  if(flag & DSPL_FLAG_FFT_SHIFT)
406  {
407  err = fft_shift(mag, n, mag);
408  if(err!=RES_OK)
409  goto error_proc;
410  }
411  }
412 
413  if(freq)
414  {
415  if(flag & DSPL_FLAG_FFT_SHIFT)
416  if(n%2)
417  err = linspace(-fs*0.5 + fs*0.5/(double)n,
418  fs*0.5 - fs*0.5/(double)n,
419  n, DSPL_SYMMETRIC, freq);
420  else
421  err = linspace(-fs*0.5, fs*0.5, n, DSPL_PERIODIC, freq);
422  else
423  err = linspace(0, fs, n, DSPL_PERIODIC, freq);
424  }
425 
426 error_proc:
427  if(X)
428  free(X);
429 
430  return err;
431 }
432 
433 
434 
435 
436 #ifdef DOXYGEN_ENGLISH
437 
438 #endif
439 #ifdef DOXYGEN_RUSSIAN
440 
441 #endif
442 int DSPL_API fft_abs_cmplx(complex_t* x, int n, fft_t* pfft,
443  double fs, int flag,
444  double* mag, double* freq)
445 {
446  int k, err = RES_OK;
447  complex_t *X = NULL;
448 
449  if(!x || !pfft)
450  return ERROR_PTR;
451 
452  if(n<1)
453  return ERROR_SIZE;
454 
455  if(mag)
456  {
457  X = (complex_t*)malloc(n*sizeof(complex_t));
458  err = fft_cmplx(x, n, pfft, X);
459  if(err!=RES_OK)
460  goto error_proc;
461 
462 
463  for(k = 0; k < n; k++)
464  mag[k] = ABS(X[k]);
465 
466  if(flag & DSPL_FLAG_FFT_SHIFT)
467  {
468  err = fft_shift(mag, n, mag);
469  if(err!=RES_OK)
470  goto error_proc;
471  }
472  }
473 
474  if(freq)
475  {
476  if(flag & DSPL_FLAG_FFT_SHIFT)
477  if(n%2)
478  err = linspace(-fs*0.5 + fs*0.5/(double)n,
479  fs*0.5 - fs*0.5/(double)n,
480  n, DSPL_SYMMETRIC, freq);
481  else
482  err = linspace(-fs*0.5, fs*0.5, n, DSPL_PERIODIC, freq);
483  else
484  err = linspace(0, fs, n, DSPL_PERIODIC, freq);
485  }
486 error_proc:
487  if(X)
488  free(X);
489 
490  return err;
491 }
492 
493 
494 
495 #ifdef DOXYGEN_ENGLISH
496 
567 #endif
568 #ifdef DOXYGEN_RUSSIAN
569 
644 #endif
645 int DSPL_API fft_cmplx(complex_t* x, int n, fft_t* pfft, complex_t* y)
646 {
647  int err;
648 
649  if(!x || !pfft || !y)
650  return ERROR_PTR;
651  if(n<1)
652  return ERROR_SIZE;
653 
654  err = fft_create(pfft, n);
655  if(err != RES_OK)
656  return err;
657 
658  memcpy(pfft->t1, x, n*sizeof(complex_t));
659 
660  return fft_krn(pfft->t1, y, pfft, n, 0);
661 }
662 
663 
664 
665 #ifdef DOXYGEN_ENGLISH
666 
667 #endif
668 #ifdef DOXYGEN_RUSSIAN
669 
670 #endif
671 int fft_krn(complex_t* t0, complex_t* t1, fft_t* p, int n, int addr)
672 {
673  int n1, n2, k, m, i;
674  complex_t *pw = p->w+addr;
675  complex_t tmp;
676 
677  n1 = 1;
678  if(n%16== 0) { n1 = 16; goto label_size; }
679  if(n%7 == 0) { n1 = 7; goto label_size; }
680  if(n%8 == 0) { n1 = 8; goto label_size; }
681  if(n%5 == 0) { n1 = 5; goto label_size; }
682  if(n%4 == 0) { n1 = 4; goto label_size; }
683  if(n%3 == 0) { n1 = 3; goto label_size; }
684  if(n%2 == 0) { n1 = 2; goto label_size; }
685 
686 label_size:
687  if(n1 == 1)
688  {
689  for(k = 0; k < n; k++)
690  {
691  RE(t1[k]) = IM(t1[k]) = 0.0;
692  for(m = 0; m < n; m++)
693  {
694  i = (k*m) % n;
695  RE(tmp) = CMRE(t0[m], pw[i]);
696  IM(tmp) = CMIM(t0[m], pw[i]);
697  RE(t1[k]) += RE(tmp);
698  IM(t1[k]) += IM(tmp);
699  }
700  }
701  }
702  else
703  {
704  n2 = n / n1;
705 
706  if(n2>1)
707  {
708  memcpy(t1, t0, n*sizeof(complex_t));
709  matrix_transpose_cmplx(t1, n2, n1, t0);
710  }
711 
712  if(n1 == 16)
713  for(k = 0; k < n2; k++)
714  dft16(t0+16*k, t1+16*k);
715 
716  if(n1 == 7)
717  for(k = 0; k < n2; k++)
718  dft7(t0+7*k, t1+7*k);
719 
720  if(n1 == 8)
721  for(k = 0; k < n2; k++)
722  dft8(t0+8*k, t1+8*k);
723 
724  if(n1 == 5)
725  for(k = 0; k < n2; k++)
726  dft5(t0+5*k, t1+5*k);
727 
728  if(n1 == 4)
729  for(k = 0; k < n2; k++)
730  dft4(t0+4*k, t1+4*k);
731 
732  if(n1 == 3)
733  for(k = 0; k < n2; k++)
734  dft3(t0+3*k, t1+3*k);
735 
736  if(n1 == 2)
737  for(k = 0; k < n2; k++)
738  dft2(t0+2*k, t1+2*k);
739 
740  if(n2 > 1)
741  {
742 
743  for(k =0; k < n; k++)
744  {
745  RE(t0[k]) = CMRE(t1[k], pw[k]);
746  IM(t0[k]) = CMIM(t1[k], pw[k]);
747  }
748 
749  matrix_transpose_cmplx(t0, n1, n2, t1);
750 
751  for(k = 0; k < n1; k++)
752  {
753  fft_krn(t1+k*n2, t0+k*n2, p, n2, addr+n);
754  }
755  matrix_transpose_cmplx(t0, n2, n1, t1);
756  }
757  }
758  return RES_OK;
759 }
760 
761 
762 
763 
764 #ifdef DOXYGEN_ENGLISH
765 
827 #endif
828 #ifdef DOXYGEN_RUSSIAN
829 
891 #endif
892 int DSPL_API fft_create(fft_t* pfft, int n)
893 {
894 
895  int n1, n2, addr, s, k, m, nw, err;
896  double phi;
897  s = n;
898  nw = addr = 0;
899 
900  if(pfft->n == n)
901  return RES_OK;
902 
903  while(s > 1)
904  {
905  n2 = 1;
906  if(s%16== 0) { n2 = 16; goto label_size; }
907  if(s%7 == 0) { n2 = 7; goto label_size; }
908  if(s%8 == 0) { n2 = 8; goto label_size; }
909  if(s%5 == 0) { n2 = 5; goto label_size; }
910  if(s%4 == 0) { n2 = 4; goto label_size; }
911  if(s%3 == 0) { n2 = 3; goto label_size; }
912  if(s%2 == 0) { n2 = 2; goto label_size; }
913 
914 
915 label_size:
916  if(n2 == 1)
917  {
918  if(s > FFT_COMPOSITE_MAX)
919  {
920  err = ERROR_FFT_SIZE;
921  goto error_proc;
922  }
923 
924  nw += s;
925  pfft->w = pfft->w ?
926  (complex_t*) realloc(pfft->w, nw*sizeof(complex_t)):
927  (complex_t*) malloc( nw*sizeof(complex_t));
928  for(k = 0; k < s; k++)
929  {
930  phi = - M_2PI * (double)k / (double)s;
931  RE(pfft->w[addr]) = cos(phi);
932  IM(pfft->w[addr]) = sin(phi);
933  addr++;
934  }
935  s = 1;
936  }
937  else
938  {
939  n1 = s / n2;
940  nw += s;
941  pfft->w = pfft->w ?
942  (complex_t*) realloc(pfft->w, nw*sizeof(complex_t)):
943  (complex_t*) malloc( nw*sizeof(complex_t));
944 
945  for(k = 0; k < n1; k++)
946  {
947  for(m = 0; m < n2; m++)
948  {
949  phi = - M_2PI * (double)(k*m) / (double)s;
950  RE(pfft->w[addr]) = cos(phi);
951  IM(pfft->w[addr]) = sin(phi);
952  addr++;
953  }
954  }
955  }
956  s /= n2;
957  }
958 
959  pfft->t0 = pfft->t0 ? (complex_t*) realloc(pfft->t0, n*sizeof(complex_t)):
960  (complex_t*) malloc( n*sizeof(complex_t));
961 
962  pfft->t1 = pfft->t1 ? (complex_t*) realloc(pfft->t1, n*sizeof(complex_t)):
963  (complex_t*) malloc( n*sizeof(complex_t));
964  pfft->n = n;
965 
966  return RES_OK;
967 error_proc:
968  if(pfft->t0) free(pfft->t0);
969  if(pfft->t1) free(pfft->t1);
970  if(pfft->w) free(pfft->w);
971  pfft->n = 0;
972  return err;
973 }
974 
975 
976 
977 
978 
979 #ifdef DOXYGEN_ENGLISH
980 
993 #endif
994 #ifdef DOXYGEN_RUSSIAN
995 
1008 #endif
1009 void DSPL_API fft_free(fft_t *pfft)
1010 {
1011  if(!pfft)
1012  return;
1013  if(pfft->w)
1014  free(pfft->w);
1015  if(pfft->t0)
1016  free(pfft->t0);
1017  if(pfft->t1)
1018  free(pfft->t1);
1019  memset(pfft, 0, sizeof(fft_t));
1020 }
1021 
1022 
1023 
1024 
1025 
1026 #ifdef DOXYGEN_ENGLISH
1027 
1028 #endif
1029 #ifdef DOXYGEN_RUSSIAN
1030 
1031 #endif
1032 int DSPL_API fft_mag(double* x, int n, fft_t* pfft,
1033  double fs, int flag,
1034  double* mag, double* freq)
1035 {
1036  int k, err;
1037  err = fft_abs(x, n, pfft, fs, flag, mag, freq);
1038  if(err != RES_OK)
1039  return err;
1040  if(mag)
1041  {
1042  if(flag & DSPL_FLAG_LOGMAG)
1043  for(k = 0; k < n; k++)
1044  mag[k] = 20.0 * log10(mag[k] + DBL_EPSILON);
1045  else
1046  for(k = 0; k < n; k++)
1047  mag[k] *= mag[k];
1048  }
1049 
1050  return err;
1051 }
1052 
1053 
1054 
1055 
1056 
1057 
1058 
1059 #ifdef DOXYGEN_ENGLISH
1060 
1061 #endif
1062 #ifdef DOXYGEN_RUSSIAN
1063 
1064 #endif
1065 int DSPL_API fft_mag_cmplx(complex_t* x, int n, fft_t* pfft,
1066  double fs, int flag,
1067  double* mag, double* freq)
1068 {
1069  int k, err;
1070  err = fft_abs_cmplx(x, n, pfft, fs, flag, mag, freq);
1071  if(err != RES_OK)
1072  return err;
1073  if(mag)
1074  {
1075  if(flag & DSPL_FLAG_LOGMAG)
1076  for(k = 0; k < n; k++)
1077  mag[k] = 20.0 * log10(mag[k] + DBL_EPSILON);
1078  else
1079  for(k = 0; k < n; k++)
1080  mag[k] *= mag[k];
1081  }
1082 
1083  return err;
1084 }
1085 
1086 
1087 
1088 
1089 
1090 #ifdef DOXYGEN_ENGLISH
1091 
1119 #endif
1120 #ifdef DOXYGEN_RUSSIAN
1121 
1152 #endif
1153 int DSPL_API fft_shift(double* x, int n, double* y)
1154 {
1155  int n2, r;
1156  int k;
1157  double tmp;
1158  double *buf;
1159 
1160  if(!x || !y)
1161  return ERROR_PTR;
1162 
1163  if(n<1)
1164  return ERROR_SIZE;
1165 
1166  r = n%2;
1167  if(!r)
1168  {
1169  n2 = n>>1;
1170  for(k = 0; k < n2; k++)
1171  {
1172  tmp = x[k];
1173  y[k] = x[k+n2];
1174  y[k+n2] = tmp;
1175  }
1176  }
1177  else
1178  {
1179  n2 = (n+1) >> 1;
1180  buf = (double*) malloc(n2*sizeof(double));
1181  memcpy(buf, x, n2*sizeof(double));
1182  memcpy(y, x+n2, (n2-1)*sizeof(double));
1183  memcpy(y+n2-1, buf, n2*sizeof(double));
1184  free(buf);
1185  }
1186  return RES_OK;
1187 }
1188 
1189 
1190 
1191 #ifdef DOXYGEN_ENGLISH
1192 
1193 #endif
1194 #ifdef DOXYGEN_RUSSIAN
1195 
1196 #endif
1197 int DSPL_API fft_shift_cmplx(complex_t* x, int n, complex_t* y)
1198 {
1199  int n2, r;
1200  int k;
1201  complex_t tmp;
1202  complex_t *buf;
1203 
1204  if(!x || !y)
1205  return ERROR_PTR;
1206 
1207  if(n<1)
1208  return ERROR_SIZE;
1209 
1210  r = n%2;
1211  if(!r)
1212  {
1213  n2 = n>>1;
1214  for(k = 0; k < n2; k++)
1215  {
1216  RE(tmp) = RE(x[k]);
1217  IM(tmp) = IM(x[k]);
1218 
1219  RE(y[k]) = RE(x[k+n2]);
1220  IM(y[k]) = IM(x[k+n2]);
1221 
1222  RE(y[k+n2]) = RE(tmp);
1223  IM(y[k+n2]) = IM(tmp);
1224  }
1225  }
1226  else
1227  {
1228  n2 = (n+1) >> 1;
1229  buf = (complex_t*) malloc(n2*sizeof(complex_t));
1230  memcpy(buf, x, n2*sizeof(complex_t));
1231  memcpy(y, x+n2, (n2-1)*sizeof(complex_t));
1232  memcpy(y+n2-1, buf, n2*sizeof(complex_t));
1233  free(buf);
1234  }
1235  return RES_OK;
1236 }
1237 
complex_t * t1
Definition: dspl.h:230
complex_t * t0
Definition: dspl.h:229
#define RE(x)
Macro sets real part of the complex number.
Definition: dspl.h:359
#define ERROR_PTR
Pointer error. This error means that one of the required pointers (memory to be allocated for) is tra...
Definition: dspl.h:549
int DSPL_API fft_shift(double *x, int n, double *y)
Perform a shift of the vector x, for use with the fft and ifft functions, in order to move the freque...
Definition: fft.c:1153
#define ERROR_SIZE
Error array size. This error occurs when in addition to the pointer the wrong input is passed to the ...
Definition: dspl.h:557
int n
Definition: dspl.h:231
int DSPL_API fft_create(fft_t *pfft, int n)
Function creates and fill fft_t structure.
Definition: fft.c:892
int DSPL_API fft_cmplx(complex_t *x, int n, fft_t *pfft, complex_t *y)
Fast Fourier transform for the complex vector.
Definition: fft.c:645
int DSPL_API ifft_cmplx(complex_t *x, int n, fft_t *pfft, complex_t *y)
Inverse fast Fourier transform.
Definition: fft.c:179
void DSPL_API fft_free(fft_t *pfft)
Free fft_t structure.
Definition: fft.c:1009
Fast Fourier Transform Object Data Structure.
Definition: dspl.h:227
int DSPL_API fft(double *x, int n, fft_t *pfft, complex_t *y)
Fast Fourier transform for the real vector.
Definition: fft.c:356
int DSPL_API re2cmplx(double *x, int n, complex_t *y)
Convert real array to the complex array.
Definition: complex.c:246
double complex_t[2]
Complex data type.
Definition: dspl.h:86
#define RES_OK
The function completed correctly. No errors.
Definition: dspl.h:497
int DSPL_API linspace(double x0, double x1, int n, int type, double *x)
Function fills a vector with n linearly spaced elements between x0 and x1.
Definition: array.c:1009
complex_t * w
Definition: dspl.h:228
#define ERROR_FFT_SIZE
The FFT size is not set correctly. FFT size can be composite , where , and – an arbitrary prime fact...
Definition: dspl.h:510
#define IM(x)
Macro sets imaginary part of the complex number.
Definition: dspl.h:417