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 <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 
378 #ifdef DOXYGEN_ENGLISH
379 
450 #endif
451 #ifdef DOXYGEN_RUSSIAN
452 
527 #endif
528 int DSPL_API fft_cmplx(complex_t* x, int n, fft_t* pfft, complex_t* y)
529 {
530  int err;
531 
532  if(!x || !pfft || !y)
533  return ERROR_PTR;
534  if(n<1)
535  return ERROR_SIZE;
536 
537  err = fft_create(pfft, n);
538  if(err != RES_OK)
539  return err;
540 
541  memcpy(pfft->t1, x, n*sizeof(complex_t));
542 
543  return fft_krn(pfft->t1, y, pfft, n, 0);
544 }
545 
546 
547 
548 #ifdef DOXYGEN_ENGLISH
549 
550 #endif
551 #ifdef DOXYGEN_RUSSIAN
552 
553 #endif
554 int fft_krn(complex_t* t0, complex_t* t1, fft_t* p, int n, int addr)
555 {
556  int n1, n2, k, m, i;
557  complex_t *pw = p->w+addr;
558  complex_t tmp;
559 
560  n1 = 1;
561  if(n%16== 0) { n1 = 16; goto label_size; }
562  if(n%7 == 0) { n1 = 7; goto label_size; }
563  if(n%8 == 0) { n1 = 8; goto label_size; }
564  if(n%5 == 0) { n1 = 5; goto label_size; }
565  if(n%4 == 0) { n1 = 4; goto label_size; }
566  if(n%3 == 0) { n1 = 3; goto label_size; }
567  if(n%2 == 0) { n1 = 2; goto label_size; }
568 
569 label_size:
570  if(n1 == 1)
571  {
572  for(k = 0; k < n; k++)
573  {
574  RE(t1[k]) = IM(t1[k]) = 0.0;
575  for(m = 0; m < n; m++)
576  {
577  i = (k*m) % n;
578  RE(tmp) = CMRE(t0[m], pw[i]);
579  IM(tmp) = CMIM(t0[m], pw[i]);
580  RE(t1[k]) += RE(tmp);
581  IM(t1[k]) += IM(tmp);
582  }
583  }
584  }
585  else
586  {
587  n2 = n / n1;
588 
589  if(n2>1)
590  {
591  memcpy(t1, t0, n*sizeof(complex_t));
592  matrix_transpose_cmplx(t1, n2, n1, t0);
593  }
594 
595  if(n1 == 16)
596  for(k = 0; k < n2; k++)
597  dft16(t0+16*k, t1+16*k);
598 
599  if(n1 == 7)
600  for(k = 0; k < n2; k++)
601  dft7(t0+7*k, t1+7*k);
602 
603  if(n1 == 8)
604  for(k = 0; k < n2; k++)
605  dft8(t0+8*k, t1+8*k);
606 
607  if(n1 == 5)
608  for(k = 0; k < n2; k++)
609  dft5(t0+5*k, t1+5*k);
610 
611  if(n1 == 4)
612  for(k = 0; k < n2; k++)
613  dft4(t0+4*k, t1+4*k);
614 
615  if(n1 == 3)
616  for(k = 0; k < n2; k++)
617  dft3(t0+3*k, t1+3*k);
618 
619  if(n1 == 2)
620  for(k = 0; k < n2; k++)
621  dft2(t0+2*k, t1+2*k);
622 
623  if(n2 > 1)
624  {
625 
626  for(k =0; k < n; k++)
627  {
628  RE(t0[k]) = CMRE(t1[k], pw[k]);
629  IM(t0[k]) = CMIM(t1[k], pw[k]);
630  }
631 
632  matrix_transpose_cmplx(t0, n1, n2, t1);
633 
634  for(k = 0; k < n1; k++)
635  {
636  fft_krn(t1+k*n2, t0+k*n2, p, n2, addr+n);
637  }
638  matrix_transpose_cmplx(t0, n2, n1, t1);
639  }
640  }
641  return RES_OK;
642 }
643 
644 
645 
646 
647 #ifdef DOXYGEN_ENGLISH
648 
710 #endif
711 #ifdef DOXYGEN_RUSSIAN
712 
774 #endif
775 int DSPL_API fft_create(fft_t* pfft, int n)
776 {
777 
778  int n1, n2, addr, s, k, m, nw, err;
779  double phi;
780  s = n;
781  nw = addr = 0;
782 
783  if(pfft->n == n)
784  return RES_OK;
785 
786  while(s > 1)
787  {
788  n2 = 1;
789  if(s%16== 0) { n2 = 16; goto label_size; }
790  if(s%7 == 0) { n2 = 7; goto label_size; }
791  if(s%8 == 0) { n2 = 8; goto label_size; }
792  if(s%5 == 0) { n2 = 5; goto label_size; }
793  if(s%4 == 0) { n2 = 4; goto label_size; }
794  if(s%3 == 0) { n2 = 3; goto label_size; }
795  if(s%2 == 0) { n2 = 2; goto label_size; }
796 
797 
798 label_size:
799  if(n2 == 1)
800  {
801  if(s > FFT_COMPOSITE_MAX)
802  {
803  err = ERROR_FFT_SIZE;
804  goto error_proc;
805  }
806 
807  nw += s;
808  pfft->w = pfft->w ?
809  (complex_t*) realloc(pfft->w, nw*sizeof(complex_t)):
810  (complex_t*) malloc( nw*sizeof(complex_t));
811  for(k = 0; k < s; k++)
812  {
813  phi = - M_2PI * (double)k / (double)s;
814  RE(pfft->w[addr]) = cos(phi);
815  IM(pfft->w[addr]) = sin(phi);
816  addr++;
817  }
818  s = 1;
819  }
820  else
821  {
822  n1 = s / n2;
823  nw += s;
824  pfft->w = pfft->w ?
825  (complex_t*) realloc(pfft->w, nw*sizeof(complex_t)):
826  (complex_t*) malloc( nw*sizeof(complex_t));
827 
828  for(k = 0; k < n1; k++)
829  {
830  for(m = 0; m < n2; m++)
831  {
832  phi = - M_2PI * (double)(k*m) / (double)s;
833  RE(pfft->w[addr]) = cos(phi);
834  IM(pfft->w[addr]) = sin(phi);
835  addr++;
836  }
837  }
838  }
839  s /= n2;
840  }
841 
842  pfft->t0 = pfft->t0 ? (complex_t*) realloc(pfft->t0, n*sizeof(complex_t)):
843  (complex_t*) malloc( n*sizeof(complex_t));
844 
845  pfft->t1 = pfft->t1 ? (complex_t*) realloc(pfft->t1, n*sizeof(complex_t)):
846  (complex_t*) malloc( n*sizeof(complex_t));
847  pfft->n = n;
848 
849  return RES_OK;
850 error_proc:
851  if(pfft->t0) free(pfft->t0);
852  if(pfft->t1) free(pfft->t1);
853  if(pfft->w) free(pfft->w);
854  pfft->n = 0;
855  return err;
856 }
857 
858 
859 
860 
861 
862 #ifdef DOXYGEN_ENGLISH
863 
876 #endif
877 #ifdef DOXYGEN_RUSSIAN
878 
891 #endif
892 void DSPL_API fft_free(fft_t *pfft)
893 {
894  if(!pfft)
895  return;
896  if(pfft->w)
897  free(pfft->w);
898  if(pfft->t0)
899  free(pfft->t0);
900  if(pfft->t1)
901  free(pfft->t1);
902  memset(pfft, 0, sizeof(fft_t));
903 }
904 
905 
906 
907 #ifdef DOXYGEN_ENGLISH
908 
909 #endif
910 #ifdef DOXYGEN_RUSSIAN
911 
912 #endif
913 int DSPL_API fft_mag(double* x, int n, fft_t* pfft,
914  double fs, int flag,
915  double* mag, double* freq)
916 {
917  int k, err = RES_OK;
918  complex_t *X = NULL;
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])+DBL_EPSILON);
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)
Макрос определяющий реальную часть комплексного числа.
Definition: dspl.h:359
#define ERROR_PTR
Ошибка указателя. Данная ошибка означает, что один из обязательных указателей (память под который дол...
Definition: dspl.h:545
int fft_shift(double *x, int n, double *y)
Перестановка спектральных отсчетов дискретного преобразования Фурье
Definition: fft.c:1098
#define ERROR_SIZE
Ошибка при передаче размера массива. Данная ошибка возникает когда помимо указателя на массив входных...
Definition: dspl.h:553
int n
Definition: dspl.h:231
int fft_create(fft_t *pfft, int n)
Заполнение структуры fft_t для алгоритма БПФ
Definition: fft.c:775
int fft_cmplx(complex_t *x, int n, fft_t *pfft, complex_t *y)
Быстрое преобразование Фурье комплексного сигнала
Definition: fft.c:528
int ifft_cmplx(complex_t *x, int n, fft_t *pfft, complex_t *y)
Обратное быстрое преобразование Фурье
Definition: fft.c:179
void fft_free(fft_t *pfft)
Очистить структуру fft_t алгоритма БПФ
Definition: fft.c:892
Структура данных объекта быстрого преобразования Фурье
Definition: dspl.h:227
int fft(double *x, int n, fft_t *pfft, complex_t *y)
Быстрое преобразование Фурье вещественного сигнала
Definition: fft.c:356
int re2cmplx(double *x, int n, complex_t *y)
Преобразование массива вещественных данных в массив комплексных данных.
Definition: complex.c:246
double complex_t[2]
Описание комплексного типа данных.
Definition: dspl.h:86
#define RES_OK
Функция завершилась корректно. Ошибки отсутствуют.
Definition: dspl.h:497
#define ABSSQR(x)
Макрос возвращает квадрат модуля комплексного числа x.
Definition: dspl.h:473
int linspace(double x0, double x1, int n, int type, double *x)
Функция заполняет массив линейно-нарастающими, равноотстоящими значениями от x0 до x1
Definition: array.c:1009
complex_t * w
Definition: dspl.h:228
#define ERROR_FFT_SIZE
Неверно задан размер БПФ. Размер БПФ может быть составным вида , где , а – произвольный простой множ...
Definition: dspl.h:510
#define IM(x)
Макрос определяющий мнимую часть комплексного числа.
Definition: dspl.h:417