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 "dspl.h"
25 #include "dspl_internal.h"
26 
27 
28 #ifdef DOXYGEN_ENGLISH
29 
102 #endif
103 #ifdef DOXYGEN_RUSSIAN
104 
177 #endif
178 int DSPL_API ifft_cmplx(complex_t *x, int n, fft_t* pfft, complex_t* y)
179 {
180  int err, k;
181  double norm;
182 
183  if(!x || !pfft || !y)
184  return ERROR_PTR;
185  if(n<1)
186  return ERROR_SIZE;
187 
188 
189  err = fft_create(pfft, n);
190  if(err != RES_OK)
191  return err;
192 
193  memcpy(pfft->t1, x, n*sizeof(complex_t));
194  for(k = 0; k < n; k++)
195  IM(pfft->t1[k]) = -IM(pfft->t1[k]);
196 
197  err = fft_krn(pfft->t1, pfft->t0, pfft, n, 0);
198 
199  if(err!=RES_OK)
200  return err;
201 
202  norm = 1.0 / (double)n;
203  for(k = 0; k < n; k++)
204  {
205  RE(y[k]) = RE(pfft->t0[k])*norm;
206  IM(y[k]) = -IM(pfft->t0[k])*norm;
207  }
208  return RES_OK;
209 }
210 
211 
212 #ifdef DOXYGEN_ENGLISH
213 
282 #endif
283 #ifdef DOXYGEN_RUSSIAN
284 
354 #endif
355 int DSPL_API fft(double* x, int n, fft_t* pfft, complex_t* y)
356 {
357  int err;
358 
359  if(!x || !pfft || !y)
360  return ERROR_PTR;
361  if(n<1)
362  return ERROR_SIZE;
363 
364 
365  err = fft_create(pfft, n);
366  if(err != RES_OK)
367  return err;
368 
369  re2cmplx(x, n, pfft->t1);
370 
371  return fft_krn(pfft->t1, y, pfft, n, 0);
372 }
373 
374 
375 
376 
377 #ifdef DOXYGEN_ENGLISH
378 
449 #endif
450 #ifdef DOXYGEN_RUSSIAN
451 
526 #endif
527 int DSPL_API fft_cmplx(complex_t* x, int n, fft_t* pfft, complex_t* y)
528 {
529  int err;
530 
531  if(!x || !pfft || !y)
532  return ERROR_PTR;
533  if(n<1)
534  return ERROR_SIZE;
535 
536  err = fft_create(pfft, n);
537  if(err != RES_OK)
538  return err;
539 
540  memcpy(pfft->t1, x, n*sizeof(complex_t));
541 
542  return fft_krn(pfft->t1, y, pfft, n, 0);
543 }
544 
545 
546 
547 #ifdef DOXYGEN_ENGLISH
548 
549 #endif
550 #ifdef DOXYGEN_RUSSIAN
551 
552 #endif
553 int fft_krn(complex_t* t0, complex_t* t1, fft_t* p, int n, int addr)
554 {
555  int n1, n2, k, m, i;
556  complex_t *pw = p->w+addr;
557  complex_t tmp;
558 
559  n1 = 1;
560  if(n%16== 0) { n1 = 16; goto label_size; }
561  if(n%7 == 0) { n1 = 7; goto label_size; }
562  if(n%8 == 0) { n1 = 8; goto label_size; }
563  if(n%5 == 0) { n1 = 5; goto label_size; }
564  if(n%4 == 0) { n1 = 4; goto label_size; }
565  if(n%3 == 0) { n1 = 3; goto label_size; }
566  if(n%2 == 0) { n1 = 2; goto label_size; }
567 
568 label_size:
569  if(n1 == 1)
570  {
571  for(k = 0; k < n; k++)
572  {
573  RE(t1[k]) = IM(t1[k]) = 0.0;
574  for(m = 0; m < n; m++)
575  {
576  i = (k*m) % n;
577  RE(tmp) = CMRE(t0[m], pw[i]);
578  IM(tmp) = CMIM(t0[m], pw[i]);
579  RE(t1[k]) += RE(tmp);
580  IM(t1[k]) += IM(tmp);
581  }
582  }
583  }
584  else
585  {
586  n2 = n / n1;
587 
588  if(n2>1)
589  {
590  memcpy(t1, t0, n*sizeof(complex_t));
591  matrix_transpose_cmplx(t1, n2, n1, t0);
592  }
593 
594  if(n1 == 16)
595  for(k = 0; k < n2; k++)
596  dft16(t0+16*k, t1+16*k);
597 
598  if(n1 == 7)
599  for(k = 0; k < n2; k++)
600  dft7(t0+7*k, t1+7*k);
601 
602  if(n1 == 8)
603  for(k = 0; k < n2; k++)
604  dft8(t0+8*k, t1+8*k);
605 
606  if(n1 == 5)
607  for(k = 0; k < n2; k++)
608  dft5(t0+5*k, t1+5*k);
609 
610  if(n1 == 4)
611  for(k = 0; k < n2; k++)
612  dft4(t0+4*k, t1+4*k);
613 
614  if(n1 == 3)
615  for(k = 0; k < n2; k++)
616  dft3(t0+3*k, t1+3*k);
617 
618  if(n1 == 2)
619  for(k = 0; k < n2; k++)
620  dft2(t0+2*k, t1+2*k);
621 
622  if(n2 > 1)
623  {
624 
625  for(k =0; k < n; k++)
626  {
627  RE(t0[k]) = CMRE(t1[k], pw[k]);
628  IM(t0[k]) = CMIM(t1[k], pw[k]);
629  }
630 
631  matrix_transpose_cmplx(t0, n1, n2, t1);
632 
633  for(k = 0; k < n1; k++)
634  {
635  fft_krn(t1+k*n2, t0+k*n2, p, n2, addr+n);
636  }
637  matrix_transpose_cmplx(t0, n2, n1, t1);
638  }
639  }
640  return RES_OK;
641 }
642 
643 
644 
645 
646 #ifdef DOXYGEN_ENGLISH
647 
709 #endif
710 #ifdef DOXYGEN_RUSSIAN
711 
773 #endif
774 int DSPL_API fft_create(fft_t* pfft, int n)
775 {
776 
777  int n1, n2, addr, s, k, m, nw, err;
778  double phi;
779  s = n;
780  nw = addr = 0;
781 
782  if(pfft->n == n)
783  return RES_OK;
784 
785  while(s > 1)
786  {
787  n2 = 1;
788  if(s%16== 0) { n2 = 16; goto label_size; }
789  if(s%7 == 0) { n2 = 7; goto label_size; }
790  if(s%8 == 0) { n2 = 8; goto label_size; }
791  if(s%5 == 0) { n2 = 5; goto label_size; }
792  if(s%4 == 0) { n2 = 4; goto label_size; }
793  if(s%3 == 0) { n2 = 3; goto label_size; }
794  if(s%2 == 0) { n2 = 2; goto label_size; }
795 
796 
797 label_size:
798  if(n2 == 1)
799  {
800  if(s > FFT_COMPOSITE_MAX)
801  {
802  err = ERROR_FFT_SIZE;
803  goto error_proc;
804  }
805 
806  nw += s;
807  pfft->w = pfft->w ?
808  (complex_t*) realloc(pfft->w, nw*sizeof(complex_t)):
809  (complex_t*) malloc( nw*sizeof(complex_t));
810  for(k = 0; k < s; k++)
811  {
812  phi = - M_2PI * (double)k / (double)s;
813  RE(pfft->w[addr]) = cos(phi);
814  IM(pfft->w[addr]) = sin(phi);
815  addr++;
816  }
817  s = 1;
818  }
819  else
820  {
821  n1 = s / n2;
822  nw += s;
823  pfft->w = pfft->w ?
824  (complex_t*) realloc(pfft->w, nw*sizeof(complex_t)):
825  (complex_t*) malloc( nw*sizeof(complex_t));
826 
827  for(k = 0; k < n1; k++)
828  {
829  for(m = 0; m < n2; m++)
830  {
831  phi = - M_2PI * (double)(k*m) / (double)s;
832  RE(pfft->w[addr]) = cos(phi);
833  IM(pfft->w[addr]) = sin(phi);
834  addr++;
835  }
836  }
837  }
838  s /= n2;
839  }
840 
841  pfft->t0 = pfft->t0 ? (complex_t*) realloc(pfft->t0, n*sizeof(complex_t)):
842  (complex_t*) malloc( n*sizeof(complex_t));
843 
844  pfft->t1 = pfft->t1 ? (complex_t*) realloc(pfft->t1, n*sizeof(complex_t)):
845  (complex_t*) malloc( n*sizeof(complex_t));
846  pfft->n = n;
847 
848  return RES_OK;
849 error_proc:
850  if(pfft->t0) free(pfft->t0);
851  if(pfft->t1) free(pfft->t1);
852  if(pfft->w) free(pfft->w);
853  pfft->n = 0;
854  return err;
855 }
856 
857 
858 
859 
860 
861 #ifdef DOXYGEN_ENGLISH
862 
875 #endif
876 #ifdef DOXYGEN_RUSSIAN
877 
890 #endif
891 void DSPL_API fft_free(fft_t *pfft)
892 {
893  if(!pfft)
894  return;
895  if(pfft->w)
896  free(pfft->w);
897  if(pfft->t0)
898  free(pfft->t0);
899  if(pfft->t1)
900  free(pfft->t1);
901  memset(pfft, 0, sizeof(fft_t));
902 }
903 
904 
905 
906 #ifdef DOXYGEN_ENGLISH
907 
908 #endif
909 #ifdef DOXYGEN_RUSSIAN
910 
911 #endif
912 int DSPL_API fft_mag(double* x, int n, fft_t* pfft,
913  double fs, int flag,
914  double* mag, double* freq)
915 {
916  int k, err = RES_OK;
917  complex_t *X = NULL;
918 
919  if(!x || !pfft)
920  return ERROR_PTR;
921 
922  if(n<1)
923  return ERROR_SIZE;
924 
925  if(mag)
926  {
927  X = (complex_t*)malloc(n*sizeof(complex_t));
928  err = fft(x, n, pfft, X);
929  if(err!=RES_OK)
930  goto error_proc;
931 
932  if(flag & DSPL_FLAG_LOGMAG)
933  for(k = 0; k < n; k++)
934  mag[k] = 10.0*log10(ABSSQR(X[k]));
935  else
936  for(k = 0; k < n; k++)
937  mag[k] = ABS(X[k]);
938  if(flag & DSPL_FLAG_FFT_SHIFT)
939  {
940  err = fft_shift(mag, n, mag);
941  if(err!=RES_OK)
942  goto error_proc;
943  }
944  }
945 
946  if(freq)
947  {
948  if(flag & DSPL_FLAG_FFT_SHIFT)
949  if(n%2)
950  err = linspace(-fs*0.5 + fs*0.5/(double)n,
951  fs*0.5 - fs*0.5/(double)n,
952  n, DSPL_SYMMETRIC, freq);
953  else
954  err = linspace(-fs*0.5, fs*0.5, n, DSPL_PERIODIC, freq);
955  else
956  err = linspace(0, fs, n, DSPL_PERIODIC, freq);
957  }
958 
959 error_proc:
960  if(X)
961  free(X);
962 
963  return err;
964 }
965 
966 
967 
968 
969 
970 
971 
972 #ifdef DOXYGEN_ENGLISH
973 
974 #endif
975 #ifdef DOXYGEN_RUSSIAN
976 
977 #endif
978 int DSPL_API fft_mag_cmplx(complex_t* x, int n, fft_t* pfft,
979  double fs, int flag,
980  double* mag, double* freq)
981 {
982  int k, err = RES_OK;
983  complex_t *X = NULL;
984 
985  if(!x || !pfft)
986  return ERROR_PTR;
987 
988  if(n<1)
989  return ERROR_SIZE;
990 
991  if(mag)
992  {
993  X = (complex_t*)malloc(n*sizeof(complex_t));
994  err = fft_cmplx(x, n, pfft, X);
995  if(err!=RES_OK)
996  goto error_proc;
997 
998  if(flag & DSPL_FLAG_LOGMAG)
999  for(k = 0; k < n; k++)
1000  mag[k] = 10.0*log10(ABSSQR(X[k]));
1001  else
1002  for(k = 0; k < n; k++)
1003  mag[k] = ABS(X[k]);
1004  if(flag & DSPL_FLAG_FFT_SHIFT)
1005  {
1006  err = fft_shift(mag, n, mag);
1007  if(err!=RES_OK)
1008  goto error_proc;
1009  }
1010  }
1011 
1012  if(freq)
1013  {
1014  if(flag & DSPL_FLAG_FFT_SHIFT)
1015  if(n%2)
1016  err = linspace(-fs*0.5 + fs*0.5/(double)n,
1017  fs*0.5 - fs*0.5/(double)n,
1018  n, DSPL_SYMMETRIC, freq);
1019  else
1020  err = linspace(-fs*0.5, fs*0.5, n, DSPL_PERIODIC, freq);
1021  else
1022  err = linspace(0, fs, n, DSPL_PERIODIC, freq);
1023  }
1024 error_proc:
1025  if(X)
1026  free(X);
1027 
1028  return err;
1029 }
1030 
1031 
1032 
1033 
1034 
1035 #ifdef DOXYGEN_ENGLISH
1036 
1064 #endif
1065 #ifdef DOXYGEN_RUSSIAN
1066 
1097 #endif
1098 int DSPL_API fft_shift(double* x, int n, double* y)
1099 {
1100  int n2, r;
1101  int k;
1102  double tmp;
1103  double *buf;
1104 
1105  if(!x || !y)
1106  return ERROR_PTR;
1107 
1108  if(n<1)
1109  return ERROR_SIZE;
1110 
1111  r = n%2;
1112  if(!r)
1113  {
1114  n2 = n>>1;
1115  for(k = 0; k < n2; k++)
1116  {
1117  tmp = x[k];
1118  y[k] = x[k+n2];
1119  y[k+n2] = tmp;
1120  }
1121  }
1122  else
1123  {
1124  n2 = (n+1) >> 1;
1125  buf = (double*) malloc(n2*sizeof(double));
1126  memcpy(buf, x, n2*sizeof(double));
1127  memcpy(y, x+n2, (n2-1)*sizeof(double));
1128  memcpy(y+n2-1, buf, n2*sizeof(double));
1129  free(buf);
1130  }
1131  return RES_OK;
1132 }
1133 
1134 
1135 
1136 #ifdef DOXYGEN_ENGLISH
1137 
1138 #endif
1139 #ifdef DOXYGEN_RUSSIAN
1140 
1141 #endif
1142 int DSPL_API fft_shift_cmplx(complex_t* x, int n, complex_t* y)
1143 {
1144  int n2, r;
1145  int k;
1146  complex_t tmp;
1147  complex_t *buf;
1148 
1149  if(!x || !y)
1150  return ERROR_PTR;
1151 
1152  if(n<1)
1153  return ERROR_SIZE;
1154 
1155  r = n%2;
1156  if(!r)
1157  {
1158  n2 = n>>1;
1159  for(k = 0; k < n2; k++)
1160  {
1161  RE(tmp) = RE(x[k]);
1162  IM(tmp) = IM(x[k]);
1163 
1164  RE(y[k]) = RE(x[k+n2]);
1165  IM(y[k]) = IM(x[k+n2]);
1166 
1167  RE(y[k+n2]) = RE(tmp);
1168  IM(y[k+n2]) = IM(tmp);
1169  }
1170  }
1171  else
1172  {
1173  n2 = (n+1) >> 1;
1174  buf = (complex_t*) malloc(n2*sizeof(complex_t));
1175  memcpy(buf, x, n2*sizeof(complex_t));
1176  memcpy(y, x+n2, (n2-1)*sizeof(complex_t));
1177  memcpy(y+n2-1, buf, n2*sizeof(complex_t));
1178  free(buf);
1179  }
1180  return RES_OK;
1181 }
1182 
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:545
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:1098
#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:553
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:774
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:527
int DSPL_API ifft_cmplx(complex_t *x, int n, fft_t *pfft, complex_t *y)
Inverse fast Fourier transform.
Definition: fft.c:178
void DSPL_API fft_free(fft_t *pfft)
Free fft_t structure.
Definition: fft.c:891
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:355
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
#define ABSSQR(x)
The macro returns the square of the modulus of a complex number x.
Definition: dspl.h:473
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