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 
26 #include "dspl.h"
27 #include "dspl_internal.h"
28 
29 
30 #ifdef DOXYGEN_ENGLISH
31 
104 #endif
105 #ifdef DOXYGEN_RUSSIAN
106 
179 #endif
180 int DSPL_API ifft_cmplx(complex_t *x, int n, fft_t* pfft, complex_t* y)
181 {
182  int err, k;
183  double norm;
184 
185  if(!x || !pfft || !y)
186  return ERROR_PTR;
187  if(n<1)
188  return ERROR_SIZE;
189 
190 
191  err = fft_create(pfft, n);
192  if(err != RES_OK)
193  return err;
194 
195  memcpy(pfft->t1, x, n*sizeof(complex_t));
196  for(k = 0; k < n; k++)
197  IM(pfft->t1[k]) = -IM(pfft->t1[k]);
198 
199  err = fft_krn(pfft->t1, pfft->t0, pfft, n, 0);
200 
201  if(err!=RES_OK)
202  return err;
203 
204  norm = 1.0 / (double)n;
205  for(k = 0; k < n; k++)
206  {
207  RE(y[k]) = RE(pfft->t0[k])*norm;
208  IM(y[k]) = -IM(pfft->t0[k])*norm;
209  }
210  return RES_OK;
211 }
212 
213 
214 #ifdef DOXYGEN_ENGLISH
215 
284 #endif
285 #ifdef DOXYGEN_RUSSIAN
286 
356 #endif
357 int DSPL_API fft(double* x, int n, fft_t* pfft, complex_t* y)
358 {
359  int err;
360 
361  if(!x || !pfft || !y)
362  return ERROR_PTR;
363  if(n<1)
364  return ERROR_SIZE;
365 
366 
367  err = fft_create(pfft, n);
368  if(err != RES_OK)
369  return err;
370 
371  re2cmplx(x, n, pfft->t1);
372 
373  return fft_krn(pfft->t1, y, pfft, n, 0);
374 }
375 
376 
377 
378 #ifdef DOXYGEN_ENGLISH
379 
380 #endif
381 #ifdef DOXYGEN_RUSSIAN
382 
383 #endif
384 int DSPL_API fft_abs(double* x, int n, fft_t* pfft,
385  double fs, int flag,
386  double* mag, double* freq)
387 {
388  int k, err = RES_OK;
389  complex_t *X = NULL;
390  if(!x || !pfft)
391  return ERROR_PTR;
392 
393  if(n<1)
394  return ERROR_SIZE;
395 
396  if(mag)
397  {
398  X = (complex_t*)malloc(n*sizeof(complex_t));
399  err = fft(x, n, pfft, X);
400  if(err!=RES_OK)
401  goto error_proc;
402 
403 
404  for(k = 0; k < n; k++)
405  mag[k] = ABS(X[k]);
406  if(flag & DSPL_FLAG_FFT_SHIFT)
407  {
408  err = fft_shift(mag, n, mag);
409  if(err!=RES_OK)
410  goto error_proc;
411  }
412  }
413 
414  if(freq)
415  {
416  if(flag & DSPL_FLAG_FFT_SHIFT)
417  if(n%2)
418  err = linspace(-fs*0.5 + fs*0.5/(double)n,
419  fs*0.5 - fs*0.5/(double)n,
420  n, DSPL_SYMMETRIC, freq);
421  else
422  err = linspace(-fs*0.5, fs*0.5, n, DSPL_PERIODIC, freq);
423  else
424  err = linspace(0, fs, n, DSPL_PERIODIC, freq);
425  }
426 
427 error_proc:
428  if(X)
429  free(X);
430 
431  return err;
432 }
433 
434 
435 
436 
437 #ifdef DOXYGEN_ENGLISH
438 
439 #endif
440 #ifdef DOXYGEN_RUSSIAN
441 
442 #endif
443 int DSPL_API fft_abs_cmplx(complex_t* x, int n, fft_t* pfft,
444  double fs, int flag,
445  double* mag, double* freq)
446 {
447  int k, err = RES_OK;
448  complex_t *X = NULL;
449 
450  if(!x || !pfft)
451  return ERROR_PTR;
452 
453  if(n<1)
454  return ERROR_SIZE;
455 
456  if(mag)
457  {
458  X = (complex_t*)malloc(n*sizeof(complex_t));
459  err = fft_cmplx(x, n, pfft, X);
460  if(err!=RES_OK)
461  goto error_proc;
462 
463 
464  for(k = 0; k < n; k++)
465  mag[k] = ABS(X[k]);
466 
467  if(flag & DSPL_FLAG_FFT_SHIFT)
468  {
469  err = fft_shift(mag, n, mag);
470  if(err!=RES_OK)
471  goto error_proc;
472  }
473  }
474 
475  if(freq)
476  {
477  if(flag & DSPL_FLAG_FFT_SHIFT)
478  if(n%2)
479  err = linspace(-fs*0.5 + fs*0.5/(double)n,
480  fs*0.5 - fs*0.5/(double)n,
481  n, DSPL_SYMMETRIC, freq);
482  else
483  err = linspace(-fs*0.5, fs*0.5, n, DSPL_PERIODIC, freq);
484  else
485  err = linspace(0, fs, n, DSPL_PERIODIC, freq);
486  }
487 error_proc:
488  if(X)
489  free(X);
490 
491  return err;
492 }
493 
494 
495 
496 #ifdef DOXYGEN_ENGLISH
497 
568 #endif
569 #ifdef DOXYGEN_RUSSIAN
570 
645 #endif
646 int DSPL_API fft_cmplx(complex_t* x, int n, fft_t* pfft, complex_t* y)
647 {
648  int err;
649 
650  if(!x || !pfft || !y)
651  return ERROR_PTR;
652  if(n<1)
653  return ERROR_SIZE;
654 
655  err = fft_create(pfft, n);
656  if(err != RES_OK)
657  return err;
658 
659  memcpy(pfft->t1, x, n*sizeof(complex_t));
660 
661  return fft_krn(pfft->t1, y, pfft, n, 0);
662 }
663 
664 
665 
666 #ifdef DOXYGEN_ENGLISH
667 
668 #endif
669 #ifdef DOXYGEN_RUSSIAN
670 
671 #endif
672 int fft_krn(complex_t* t0, complex_t* t1, fft_t* p, int n, int addr)
673 {
674  int n1, n2, k, m, i;
675  complex_t *pw = p->w+addr;
676  complex_t tmp;
677 
678  n1 = 1;
679  if(n % 4096 == 0) { n1 = 4096; goto label_size; }
680  if(n % 2048 == 0) { n1 = 2048; goto label_size; }
681  if(n % 1024 == 0) { n1 = 1024; goto label_size; }
682  if(n % 512 == 0) { n1 = 512; goto label_size; }
683  if(n % 256 == 0) { n1 = 256; goto label_size; }
684  if(n % 128 == 0) { n1 = 128; goto label_size; }
685  if(n % 64 == 0) { n1 = 64; goto label_size; }
686  if(n % 32 == 0) { n1 = 32; goto label_size; }
687  if(n % 16 == 0) { n1 = 16; goto label_size; }
688  if(n % 7 == 0) { n1 = 7; goto label_size; }
689  if(n % 8 == 0) { n1 = 8; goto label_size; }
690  if(n % 5 == 0) { n1 = 5; goto label_size; }
691  if(n % 4 == 0) { n1 = 4; goto label_size; }
692  if(n % 3 == 0) { n1 = 3; goto label_size; }
693  if(n % 2 == 0) { n1 = 2; goto label_size; }
694 
695 label_size:
696  if(n1 == 1)
697  {
698  for(k = 0; k < n; k++)
699  {
700  RE(t1[k]) = IM(t1[k]) = 0.0;
701  for(m = 0; m < n; m++)
702  {
703  i = (k*m) % n;
704  RE(tmp) = CMRE(t0[m], pw[i]);
705  IM(tmp) = CMIM(t0[m], pw[i]);
706  RE(t1[k]) += RE(tmp);
707  IM(t1[k]) += IM(tmp);
708  }
709  }
710  }
711  else
712  {
713  n2 = n / n1;
714 
715  if(n2>1)
716  {
717  memcpy(t1, t0, n*sizeof(complex_t));
718  matrix_transpose_cmplx(t1, n2, n1, t0);
719  }
720 
721  if(n1 == 4096)
722  for(k = 0; k < n2; k++)
723  dft4096(t0+4096*k, t1+4096*k, p->w4096, p->w256);
724 
725  if(n1 == 2048)
726  for(k = 0; k < n2; k++)
727  dft2048(t0+2048*k, t1+2048*k, p->w2048, p->w32, p->w64);
728 
729  if(n1 == 1024)
730  for(k = 0; k < n2; k++)
731  dft1024(t0+1024*k, t1+1024*k, p->w1024, p->w32);
732 
733  if(n1 == 512)
734  for(k = 0; k < n2; k++)
735  dft512(t0+512*k, t1+512*k, p->w512, p->w32);
736 
737  if(n1 == 256)
738  for(k = 0; k < n2; k++)
739  dft256(t0+256*k, t1+256*k, p->w256);
740 
741  if(n1 == 128)
742  for(k = 0; k < n2; k++)
743  dft128(t0+128*k, t1+128*k, p->w128);
744 
745  if(n1 == 64)
746  for(k = 0; k < n2; k++)
747  dft64(t0+64*k, t1+64*k, p->w64);
748 
749  if(n1 == 32)
750  for(k = 0; k < n2; k++)
751  dft32(t0+32*k, t1+32*k, p->w32);
752 
753  if(n1 == 16)
754  for(k = 0; k < n2; k++)
755  dft16(t0+16*k, t1+16*k);
756 
757  if(n1 == 7)
758  for(k = 0; k < n2; k++)
759  dft7(t0+7*k, t1+7*k);
760 
761  if(n1 == 8)
762  for(k = 0; k < n2; k++)
763  dft8(t0+8*k, t1+8*k);
764 
765  if(n1 == 5)
766  for(k = 0; k < n2; k++)
767  dft5(t0+5*k, t1+5*k);
768 
769  if(n1 == 4)
770  for(k = 0; k < n2; k++)
771  dft4(t0+4*k, t1+4*k);
772 
773  if(n1 == 3)
774  for(k = 0; k < n2; k++)
775  dft3(t0+3*k, t1+3*k);
776 
777  if(n1 == 2)
778  for(k = 0; k < n2; k++)
779  dft2(t0+2*k, t1+2*k);
780 
781  if(n2 > 1)
782  {
783 
784  for(k =0; k < n; k++)
785  {
786  RE(t0[k]) = CMRE(t1[k], pw[k]);
787  IM(t0[k]) = CMIM(t1[k], pw[k]);
788  }
789 
790  matrix_transpose_cmplx(t0, n1, n2, t1);
791 
792 
793  for(k = 0; k < n1; k++)
794  {
795  fft_krn(t1+k*n2, t0+k*n2, p, n2, addr+n);
796  }
797 
798  matrix_transpose_cmplx(t0, n2, n1, t1);
799  }
800  }
801  return RES_OK;
802 }
803 
804 
805 
806 
807 #ifdef DOXYGEN_ENGLISH
808 
870 #endif
871 #ifdef DOXYGEN_RUSSIAN
872 
934 #endif
935 int DSPL_API fft_create(fft_t* pfft, int n)
936 {
937 
938  int n1, n2, addr, s, k, m, nw, err;
939  double phi;
940  s = n;
941  nw = addr = 0;
942 
943  if(pfft->n == n)
944  return RES_OK;
945 
946  while(s > 1)
947  {
948  n2 = 1;
949  if(s%4096 == 0) { n2 = 4096; goto label_size; }
950  if(s%2048 == 0) { n2 = 2048; goto label_size; }
951  if(s%1024 == 0) { n2 = 1024; goto label_size; }
952  if(s%512 == 0) { n2 = 512; goto label_size; }
953  if(s%256 == 0) { n2 = 256; goto label_size; }
954  if(s%128 == 0) { n2 = 128; goto label_size; }
955  if(s% 64 == 0) { n2 = 64; goto label_size; }
956  if(s% 32 == 0) { n2 = 32; goto label_size; }
957  if(s% 16 == 0) { n2 = 16; goto label_size; }
958  if(s% 7 == 0) { n2 = 7; goto label_size; }
959  if(s% 8 == 0) { n2 = 8; goto label_size; }
960  if(s% 5 == 0) { n2 = 5; goto label_size; }
961  if(s% 4 == 0) { n2 = 4; goto label_size; }
962  if(s% 3 == 0) { n2 = 3; goto label_size; }
963  if(s% 2 == 0) { n2 = 2; goto label_size; }
964 
965 
966 label_size:
967  if(n2 == 1)
968  {
969  if(s > FFT_COMPOSITE_MAX)
970  {
971  err = ERROR_FFT_SIZE;
972  goto error_proc;
973  }
974 
975  nw += s;
976  pfft->w = pfft->w ?
977  (complex_t*) realloc(pfft->w, nw*sizeof(complex_t)):
978  (complex_t*) malloc( nw*sizeof(complex_t));
979  for(k = 0; k < s; k++)
980  {
981  phi = - M_2PI * (double)k / (double)s;
982  RE(pfft->w[addr]) = cos(phi);
983  IM(pfft->w[addr]) = sin(phi);
984  addr++;
985  }
986  s = 1;
987  }
988  else
989  {
990  n1 = s / n2;
991  nw += s;
992  pfft->w = pfft->w ?
993  (complex_t*) realloc(pfft->w, nw*sizeof(complex_t)):
994  (complex_t*) malloc( nw*sizeof(complex_t));
995 
996  for(k = 0; k < n1; k++)
997  {
998  for(m = 0; m < n2; m++)
999  {
1000  phi = - M_2PI * (double)(k*m) / (double)s;
1001  RE(pfft->w[addr]) = cos(phi);
1002  IM(pfft->w[addr]) = sin(phi);
1003  addr++;
1004  }
1005  }
1006  }
1007  s /= n2;
1008  }
1009 
1010  pfft->t0 = pfft->t0 ? (complex_t*) realloc(pfft->t0, n*sizeof(complex_t)):
1011  (complex_t*) malloc( n*sizeof(complex_t));
1012 
1013  pfft->t1 = pfft->t1 ? (complex_t*) realloc(pfft->t1, n*sizeof(complex_t)):
1014  (complex_t*) malloc( n*sizeof(complex_t));
1015  pfft->n = n;
1016 
1017  /* w32 fill */
1018  addr = 0;
1019  for(k = 0; k < 4; k++)
1020  {
1021  for(m = 0; m < 8; m++)
1022  {
1023  phi = - M_2PI * (double)(k*m) / 32.0;
1024  RE(pfft->w32[addr]) = cos(phi);
1025  IM(pfft->w32[addr]) = sin(phi);
1026  addr++;
1027  }
1028  }
1029 
1030 
1031  /* w64 fill */
1032  addr = 0;
1033  for(k = 0; k < 8; k++)
1034  {
1035  for(m = 0; m < 8; m++)
1036  {
1037  phi = - M_2PI * (double)(k*m) / 64.0;
1038  RE(pfft->w64[addr]) = cos(phi);
1039  IM(pfft->w64[addr]) = sin(phi);
1040  addr++;
1041  }
1042  }
1043 
1044  /* w128 fill */
1045  addr = 0;
1046  for(k = 0; k < 8; k++)
1047  {
1048  for(m = 0; m < 16; m++)
1049  {
1050  phi = - M_2PI * (double)(k*m) / 128.0;
1051  RE(pfft->w128[addr]) = cos(phi);
1052  IM(pfft->w128[addr]) = sin(phi);
1053  addr++;
1054  }
1055  }
1056 
1057  /* w256 fill */
1058  addr = 0;
1059  for(k = 0; k < 16; k++)
1060  {
1061  for(m = 0; m < 16; m++)
1062  {
1063  phi = - M_2PI * (double)(k*m) / 256.0;
1064  RE(pfft->w256[addr]) = cos(phi);
1065  IM(pfft->w256[addr]) = sin(phi);
1066  addr++;
1067  }
1068  }
1069 
1070  /* w512 fill */
1071  addr = 0;
1072  for(k = 0; k < 16; k++)
1073  {
1074  for(m = 0; m < 32; m++)
1075  {
1076  phi = - M_2PI * (double)(k*m) / 512.0;
1077  RE(pfft->w512[addr]) = cos(phi);
1078  IM(pfft->w512[addr]) = sin(phi);
1079  addr++;
1080  }
1081  }
1082 
1083  /* w1024 fill */
1084  if(pfft->w1024 == NULL)
1085  {
1086  pfft->w1024 = (complex_t*) malloc(1024 * sizeof(complex_t));
1087  addr = 0;
1088  for(k = 0; k < 32; k++)
1089  {
1090  for(m = 0; m < 32; m++)
1091  {
1092  phi = - M_2PI * (double)(k*m) / 1024.0;
1093  RE(pfft->w1024[addr]) = cos(phi);
1094  IM(pfft->w1024[addr]) = sin(phi);
1095  addr++;
1096  }
1097  }
1098  }
1099 
1100  /* w2048 fill */
1101  if(pfft->w2048 == NULL)
1102  {
1103  pfft->w2048= (complex_t*) malloc(2048 * sizeof(complex_t));
1104  addr = 0;
1105  for(k = 0; k < 32; k++)
1106  {
1107  for(m = 0; m < 64; m++)
1108  {
1109  phi = - M_2PI * (double)(k*m) / 2048.0;
1110  RE(pfft->w2048[addr]) = cos(phi);
1111  IM(pfft->w2048[addr]) = sin(phi);
1112  addr++;
1113  }
1114  }
1115  }
1116 
1117  /* w4096 fill */
1118  if(pfft->w4096 == NULL)
1119  {
1120  pfft->w4096= (complex_t*) malloc(4096 * sizeof(complex_t));
1121  addr = 0;
1122  for(k = 0; k < 16; k++)
1123  {
1124  for(m = 0; m < 256; m++)
1125  {
1126  phi = - M_2PI * (double)(k*m) / 4096.0;
1127  RE(pfft->w4096[addr]) = cos(phi);
1128  IM(pfft->w4096[addr]) = sin(phi);
1129  addr++;
1130  }
1131  }
1132  }
1133 
1134  return RES_OK;
1135 error_proc:
1136  if(pfft->t0) free(pfft->t0);
1137  if(pfft->t1) free(pfft->t1);
1138  if(pfft->w) free(pfft->w);
1139  pfft->n = 0;
1140  return err;
1141 }
1142 
1143 
1144 
1145 
1146 
1147 #ifdef DOXYGEN_ENGLISH
1148 
1161 #endif
1162 #ifdef DOXYGEN_RUSSIAN
1163 
1176 #endif
1177 void DSPL_API fft_free(fft_t *pfft)
1178 {
1179  if(!pfft)
1180  return;
1181  if(pfft->w)
1182  free(pfft->w);
1183  if(pfft->t0)
1184  free(pfft->t0);
1185  if(pfft->t1)
1186  free(pfft->t1);
1187 
1188  if(pfft->w1024)
1189  free(pfft->w1024);
1190 
1191  if(pfft->w2048)
1192  free(pfft->w2048);
1193 
1194  if(pfft->w4096)
1195  free(pfft->w4096);
1196 
1197  memset(pfft, 0, sizeof(fft_t));
1198 }
1199 
1200 
1201 
1202 
1203 
1204 #ifdef DOXYGEN_ENGLISH
1205 
1206 #endif
1207 #ifdef DOXYGEN_RUSSIAN
1208 
1209 #endif
1210 int DSPL_API fft_mag(double* x, int n, fft_t* pfft,
1211  double fs, int flag,
1212  double* mag, double* freq)
1213 {
1214  int k, err;
1215  fft_t *cfft = NULL;
1216 
1217  if(pfft)
1218  cfft = pfft;
1219  else
1220  {
1221  cfft = (fft_t*) malloc(sizeof(fft_t));
1222  memset(cfft, 0, sizeof(fft_t));
1223  }
1224 
1225  err = fft_abs(x, n, cfft, fs, flag, mag, freq);
1226  if(err != RES_OK)
1227  goto error_proc;
1228 
1229  if(mag)
1230  {
1231  if(flag & DSPL_FLAG_LOGMAG)
1232  for(k = 0; k < n; k++)
1233  mag[k] = 20.0 * log10(mag[k] + DBL_EPSILON);
1234  else
1235  for(k = 0; k < n; k++)
1236  mag[k] *= mag[k];
1237  }
1238 error_proc:
1239  if(cfft && cfft != pfft)
1240  free(cfft);
1241  return err;
1242 }
1243 
1244 
1245 
1246 
1247 
1248 
1249 
1250 #ifdef DOXYGEN_ENGLISH
1251 
1252 #endif
1253 #ifdef DOXYGEN_RUSSIAN
1254 
1255 #endif
1256 int DSPL_API fft_mag_cmplx(complex_t* x, int n, fft_t* pfft,
1257  double fs, int flag,
1258  double* mag, double* freq)
1259 {
1260  int k, err;
1261  fft_t *cfft = NULL;
1262 
1263  if(pfft)
1264  cfft = pfft;
1265  else
1266  {
1267  cfft = (fft_t*) malloc(sizeof(fft_t));
1268  memset(cfft, 0, sizeof(fft_t));
1269  }
1270 
1271  err = fft_abs_cmplx(x, n, cfft, fs, flag, mag, freq);
1272  if(err != RES_OK)
1273  goto error_proc;
1274 
1275  if(mag)
1276  {
1277  if(flag & DSPL_FLAG_LOGMAG)
1278  for(k = 0; k < n; k++)
1279  mag[k] = 20.0 * log10(mag[k] + DBL_EPSILON);
1280  else
1281  for(k = 0; k < n; k++)
1282  mag[k] *= mag[k];
1283  }
1284 error_proc:
1285  if(cfft && cfft != pfft)
1286  free(cfft);
1287  return err;
1288 }
1289 
1290 
1291 
1292 
1293 
1294 #ifdef DOXYGEN_ENGLISH
1295 
1323 #endif
1324 #ifdef DOXYGEN_RUSSIAN
1325 
1356 #endif
1357 int DSPL_API fft_shift(double* x, int n, double* y)
1358 {
1359  int n2, r;
1360  int k;
1361  double tmp;
1362  double *buf;
1363 
1364  if(!x || !y)
1365  return ERROR_PTR;
1366 
1367  if(n<1)
1368  return ERROR_SIZE;
1369 
1370  r = n%2;
1371  if(!r)
1372  {
1373  n2 = n>>1;
1374  for(k = 0; k < n2; k++)
1375  {
1376  tmp = x[k];
1377  y[k] = x[k+n2];
1378  y[k+n2] = tmp;
1379  }
1380  }
1381  else
1382  {
1383  n2 = (n+1) >> 1;
1384  buf = (double*) malloc(n2*sizeof(double));
1385  memcpy(buf, x, n2*sizeof(double));
1386  memcpy(y, x+n2, (n2-1)*sizeof(double));
1387  memcpy(y+n2-1, buf, n2*sizeof(double));
1388  free(buf);
1389  }
1390  return RES_OK;
1391 }
1392 
1393 
1394 
1395 #ifdef DOXYGEN_ENGLISH
1396 
1397 #endif
1398 #ifdef DOXYGEN_RUSSIAN
1399 
1400 #endif
1401 int DSPL_API fft_shift_cmplx(complex_t* x, int n, complex_t* y)
1402 {
1403  int n2, r;
1404  int k;
1405  complex_t tmp;
1406  complex_t *buf;
1407 
1408  if(!x || !y)
1409  return ERROR_PTR;
1410 
1411  if(n<1)
1412  return ERROR_SIZE;
1413 
1414  r = n%2;
1415  if(!r)
1416  {
1417  n2 = n>>1;
1418  for(k = 0; k < n2; k++)
1419  {
1420  RE(tmp) = RE(x[k]);
1421  IM(tmp) = IM(x[k]);
1422 
1423  RE(y[k]) = RE(x[k+n2]);
1424  IM(y[k]) = IM(x[k+n2]);
1425 
1426  RE(y[k+n2]) = RE(tmp);
1427  IM(y[k+n2]) = IM(tmp);
1428  }
1429  }
1430  else
1431  {
1432  n2 = (n+1) >> 1;
1433  buf = (complex_t*) malloc(n2*sizeof(complex_t));
1434  memcpy(buf, x, n2*sizeof(complex_t));
1435  memcpy(y, x+n2, (n2-1)*sizeof(complex_t));
1436  memcpy(y+n2-1, buf, n2*sizeof(complex_t));
1437  free(buf);
1438  }
1439  return RES_OK;
1440 }
1441 
complex_t * w1024
Definition: dspl.h:289
complex_t * t1
Definition: dspl.h:281
complex_t * t0
Definition: dspl.h:280
#define RE(x)
Macro sets real part of the complex number.
Definition: dspl.h:420
#define ERROR_PTR
Pointer error. This error means that one of the required pointers (memory to be allocated for) is tra...
Definition: dspl.h:610
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:1357
#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:618
int n
Definition: dspl.h:292
complex_t w512[512]
Definition: dspl.h:288
int DSPL_API fft_create(fft_t *pfft, int n)
Function creates and fill fft_t structure.
Definition: fft.c:935
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:646
int DSPL_API ifft_cmplx(complex_t *x, int n, fft_t *pfft, complex_t *y)
Inverse fast Fourier transform.
Definition: fft.c:180
void DSPL_API fft_free(fft_t *pfft)
Free fft_t structure.
Definition: fft.c:1177
Fast Fourier Transform Object Data Structure.
Definition: dspl.h:278
complex_t * w4096
Definition: dspl.h:291
int DSPL_API fft(double *x, int n, fft_t *pfft, complex_t *y)
Fast Fourier transform for the real vector.
Definition: fft.c:357
int DSPL_API re2cmplx(double *x, int n, complex_t *y)
Convert real array to the complex array.
Definition: complex.c:246
complex_t w64[64]
Definition: dspl.h:285
double complex_t[2]
Complex data type.
Definition: dspl.h:86
complex_t w32[32]
Definition: dspl.h:284
#define RES_OK
The function completed correctly. No errors.
Definition: dspl.h:558
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:1041
complex_t * w
Definition: dspl.h:279
#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:571
complex_t * w2048
Definition: dspl.h:290
complex_t w128[128]
Definition: dspl.h:286
complex_t w256[256]
Definition: dspl.h:287
#define IM(x)
Macro sets imaginary part of the complex number.
Definition: dspl.h:478