added verificator and FFT7. Also added readbin

Makefile

@@ -24,6 +24,8 @@ all:
 	cp  -r include/dspl.c release/include/dspl.c
 	$(MAKE) -f Makefile.test
 	cp  -r $(RELEASE_DIR)/$(DSPL_LIBNAME) test/bin/$(DSPL_LIBNAME)
+	$(MAKE) -f Makefile.verif
+	cp  -r $(RELEASE_DIR)/$(DSPL_LIBNAME) verif/bin/$(DSPL_LIBNAME)
 
 	
 	

Makefile.verif

@@ -0,0 +1,28 @@
+PRJ_DIR = verif
+SRC_DIR = $(PRJ_DIR)/src
+BIN_DIR = $(PRJ_DIR)/bin
+
+
+include Makefile.dirs
+
+DSPL_C_FILE = $(INC_DIR)/dspl.c
+DSPL_O_FILE = $(PRJ_DIR)/obj/dspl.o
+
+SRC_FILES	= $(wildcard $(SRC_DIR)/*.c)
+BIN_FILES	= $(addprefix $(BIN_DIR)/,$(notdir $(SRC_FILES:.c=.exe)))
+
+CFLAGS  = -Wall -O3 -I$(INC_DIR) -D$(DEF_OS) 
+
+all: $(BIN_FILES)
+
+$(BIN_DIR)/%.exe:$(SRC_DIR)/%.c $(DSPL_O_FILE)
+	$(CC)  $(CFLAGS)  $< $(DSPL_O_FILE)  -o $@ $(LFLAGS) 
+
+$(DSPL_O_FILE):$(DSPL_C_FILE)
+	$(CC) -c $(CFLAGS) $(DSPL_C_FILE) -o  $(DSPL_O_FILE) $(LFLAGS)
+
+clean:
+	rm -f   $(BIN_DIR)/*.exe
+	rm -f   $(BIN_DIR)/$(DSPL_LIBNAME)
+	rm -f	$(DSPL_O_FILE)
+

dspl/src/array.c

@@ -152,3 +152,87 @@ int DSPL_API flipip_cmplx(complex_t* x, int n)
   return RES_OK;
 }
 
+
+
+
+/******************************************************************************
+Verif double
+*******************************************************************************/
+int DSPL_API verif(double* x,  double* y, size_t n, double eps, double* err)
+{
+  double d, maxd; 
+  size_t k; 
+  if(!x || !y)
+    return ERROR_PTR;
+  if(n < 1)
+    return ERROR_SIZE;
+  if(eps <= 0.0 )
+    return ERROR_NEGATIVE;
+    
+  maxd = -100.0;
+  
+  for(k = 0; k < n; k++)
+  {
+    d = fabs(x[k] - y[k]);
+    if(fabs(x[k]) > 0.0)
+    {
+      d = d / fabs(x[k]);
+      if(d > maxd)
+        maxd = d;
+    }
+  }
+  if(err) 
+    *err = maxd;
+    
+  if(maxd > eps)
+    err = DSPL_VERIF_FAILED;
+  else
+    err = DSPL_VERIF_SUCCESS;
+ 
+  return err;
+}
+
+
+
+/******************************************************************************
+Verif double
+*******************************************************************************/
+int DSPL_API verif_cmplx(complex_t* x,  complex_t* y, size_t n, 
+                         double eps, double* err)
+{
+  
+  complex_t d;
+  double mx, md, maxd; 
+  size_t k; 
+  if(!x || !y)
+    return ERROR_PTR;
+  if(n < 1)
+    return ERROR_SIZE;
+  if(eps <= 0.0 )
+    return ERROR_NEGATIVE;
+    
+  maxd = -100.0;
+  
+  for(k = 0; k < n; k++)
+  {
+    RE(d) = RE(x[k]) - RE(y[k]);
+    IM(d) = IM(x[k]) - IM(y[k]);
+    md = ABS(d);
+    mx = ABS(x[k]);
+    if(mx > 0.0)
+    {
+      md = md / mx;
+      if(md > maxd)
+        maxd = md;
+    }
+  }
+  if(err) 
+    *err = maxd;
+    
+  if(maxd > eps)
+    err = DSPL_VERIF_FAILED;
+  else
+    err = DSPL_VERIF_SUCCESS;
+ 
+  return err;
+}

dspl/src/dspl_internal.h

@@ -18,7 +18,40 @@ void fft_dit_krn(complex_t *x0, complex_t *x1, complex_t *w, int n,
                  complex_t *y0, complex_t *y1);
 
 int fft_p2(int n);
+
+
 void dft2 (complex_t *x,  complex_t* y);
 
+#define DFT3_W         0.866025403784439  
+void dft3 (complex_t *x,  complex_t* y);
+
+void dft4 (complex_t *x,  complex_t* y);
+
+#define DFT5_W1       -1.250000000000000         
+#define DFT5_W2        0.559016994374947         
+#define DFT5_W3        1.538841768587630         
+#define DFT5_W4        0.587785252292473        
+#define DFT5_W5        0.363271264002680      
+void dft5 (complex_t *x,  complex_t* y);
+
+#define DFT7_W1      -1.166666666666666518636930
+#define DFT7_W2       0.790156468525400224045541
+#define DFT7_W3       0.055854267289647742400494
+#define DFT7_W4       0.734302201235752405317214
+#define DFT7_W5      -0.440958551844098378680314
+#define DFT7_W6      -0.340872930623931369442658
+#define DFT7_W7       0.533969360337725240661655
+#define DFT7_W8      -0.874842290961656665615465    
+void dft7 (complex_t *x,  complex_t* y);
+
+#define DFT16_W1       0.923879532511287
+#define DFT16_W2       0.382683432365090
+#define DFT16_W3       0.707106781186548
+void dft16 (complex_t *x,  complex_t* y);
+
+void transpose16x16(complex_t *x,  complex_t* y);
+
+
+
 #endif
 

dspl/src/fft.c

@@ -24,25 +24,6 @@
 #include "dspl.h"
 #include "dspl_internal.h"
 
-
-
-
-
-/*******************************************************************************
-2 points DFT
-*******************************************************************************/
-void dft2 (complex_t *x,  complex_t* y)
-{
-  RE(y[0]) = RE(x[0]) + RE(x[1]);
-  IM(y[0]) = IM(x[0]) + IM(x[1]);
-
-  RE(y[1]) = RE(x[0]) - RE(x[1]);
-  IM(y[1]) = IM(x[0]) - IM(x[1]);
-}
-
-
-
-
 /*******************************************************************************
 COMPLEX vector IFFT
 *******************************************************************************/
@@ -80,11 +61,6 @@ int DSPL_API ifft_cmplx(complex_t *x, int n, fft_t* pfft, complex_t* y)
 }
 
 
-
-
-
-
-
 /*******************************************************************************
 Real vector FFT
 *******************************************************************************/
@@ -257,29 +233,50 @@ int DSPL_API fftn_krn(complex_t* t0, complex_t* t1, fft_t* p, int n, int addr)
   complex_t *pw = p->w+addr;
 
   n1 = 1;
-  if(n%2 == 0)
-  {
-    n1 = 2;
-    goto label_size;
-  }
+  if(n%16== 0) { n1 = 16;   goto label_size;  }
+  if(n%7 == 0) { n1 =  7;   goto label_size;  }
+  if(n%5 == 0) { n1 =  5;   goto label_size;  }
+  if(n%4 == 0) { n1 =  4;   goto label_size;  }
+  if(n%3 == 0) { n1 =  3;   goto label_size;  }
+  if(n%2 == 0) { n1 =  2;   goto label_size;  }
 
 label_size:
   if(n1 == 1)
     return ERROR_FFT_SIZE;
+  
   n2 = n / n1;
-  memcpy(t1, t0, n*sizeof(complex_t));
-
-  transpose_cmplx(t1, n2, n1, t0);
-
-  if(n1 == 2)
+  
+  if(n2>1)
   {
-    for(k = 0; k < n2; k++)
-    {
-      dft2(t0+2*k, t1+2*k);
-    }
+    memcpy(t1, t0, n*sizeof(complex_t));
+   transpose_cmplx(t1, n2, n1, t0);
   }
+ 
+  if(n1 == 16)
+    for(k = 0; k < n2; k++)
+      dft16(t0+16*k, t1+16*k);
+      
+  if(n1 == 7)
+    for(k = 0; k < n2; k++)
+      dft7(t0+7*k, t1+7*k); 
+      
+  if(n1 == 5)
+    for(k = 0; k < n2; k++)
+      dft5(t0+5*k, t1+5*k);
+ 
+  if(n1 == 4)
+    for(k = 0; k < n2; k++)
+      dft4(t0+4*k, t1+4*k);
+  
+  if(n1 == 3)
+    for(k = 0; k < n2; k++)
+      dft3(t0+3*k, t1+3*k);
+  
+  if(n1 == 2)
+    for(k = 0; k < n2; k++)
+      dft2(t0+2*k, t1+2*k);
 
-  if(n1 > 1)
+  if(n2 > 1)
   {
 
     for(k =0; k < n; k++)
@@ -289,14 +286,11 @@ label_size:
     }
 
     transpose_cmplx(t0, n1, n2, t1);
-
-
+    
     for(k = 0; k < n1; k++)
     {
       fftn_krn(t1+k*n2, t0+k*n2, p, n2, addr+n);
     }
-
-
     transpose_cmplx(t0, n2, n1, t1);
   }
   return RES_OK;
@@ -320,16 +314,18 @@ int DSPL_API fftn_create(fft_t *pfft, int n)
   while(s > 1)
   {
     n2 = 1;
-    if(s%2 == 0)
-    {
-      n2 = 2;
-      goto label_size;
-    }
+    if(s%16== 0) { n2 = 16; goto label_size; }
+    if(s%7 == 0) { n2 =  7; goto label_size; }
+    if(s%5 == 0) { n2 =  5; goto label_size; }
+    if(s%4 == 0) { n2 =  4; goto label_size; }
+    if(s%3 == 0) { n2 =  3; goto label_size; }
+    if(s%2 == 0) { n2 =  2; goto label_size; }
 
 
 label_size:
     if(n2 == 1)
       return ERROR_FFT_SIZE;
+    
     n1 = s / n2;
     nw += s;
     pfft->w = pfft->w ? (complex_t*) realloc(pfft->w,  nw*sizeof(complex_t)):

dspl/src/fft_subkernel.c

@@ -0,0 +1,417 @@
+/*
+* Copyright (c) 2015-2018 Sergey Bakhurin
+* Digital Signal Processing Library [http://dsplib.org]
+*
+* This file is part of libdspl-2.0.
+*
+* is free software: you can redistribute it and/or modify
+* it under the terms of the GNU Lesser  General Public License as published by
+* the Free Software Foundation, either version 3 of the License, or
+* (at your option) any later version.
+*
+* DSPL is distributed in the hope that it will be useful,
+* but WITHOUT ANY WARRANTY; without even the implied warranty of
+* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+* GNU General Public License for more details.
+*
+* You should have received a copy of the GNU Lesser General Public License
+* along with Foobar.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#include <stdlib.h>
+#include <stdio.h>
+#include <string.h>
+#include "dspl.h"
+#include "dspl_internal.h"
+
+
+
+/*******************************************************************************
+2 points DFT
+*******************************************************************************/
+void dft2 (complex_t *x,  complex_t* y)
+{
+  RE(y[0]) = RE(x[0]) + RE(x[1]);
+  IM(y[0]) = IM(x[0]) + IM(x[1]);
+
+  RE(y[1]) = RE(x[0]) - RE(x[1]);
+  IM(y[1]) = IM(x[0]) - IM(x[1]);
+}
+
+
+
+/*******************************************************************************
+3 points DFT (Winograd algorithm)
+*******************************************************************************/
+void dft3 (complex_t *x,  complex_t* y)
+{
+  double a, b, c, d;
+  
+  a = RE(x[0]) - 0.5 * (RE(x[1]) + RE(x[2]));
+  b =         DFT3_W * (IM(x[1]) - IM(x[2]));
+  c = IM(x[0]) - 0.5 * (IM(x[1]) + IM(x[2]));
+  d =         DFT3_W * (RE(x[2]) - RE(x[1]));
+   
+  RE(y[0]) = RE(x[0]) + RE(x[1]) + RE(x[2]);
+  IM(y[0]) = IM(x[0]) + IM(x[1]) + IM(x[2]);
+  
+  RE(y[1]) = a + b;
+  IM(y[1]) = c + d;
+
+  RE(y[2]) = a - b;
+  IM(y[2]) = c - d;
+}
+
+
+
+/*******************************************************************************
+4 points DFT (Winograd algorithm)
+*******************************************************************************/
+void dft4 (complex_t *x,  complex_t* y)
+{
+  double a0, b0, c0, d0, a1, b1, c1, d1;
+  
+  a0 = RE(x[0]) + RE(x[2]);
+  b0 = RE(x[1]) + RE(x[3]);
+  c0 = IM(x[0]) + IM(x[2]);
+  d0 = IM(x[1]) + IM(x[3]);
+  
+  a1 = RE(x[0]) - RE(x[2]);
+  b1 = RE(x[1]) - RE(x[3]);
+  c1 = IM(x[0]) - IM(x[2]);
+  d1 = IM(x[1]) - IM(x[3]);
+  
+  RE(y[0]) = a0 + b0;
+  IM(y[0]) = c0 + d0;
+
+  RE(y[1]) = a1 + d1;
+  IM(y[1]) = c1 - b1;
+  
+  RE(y[2]) = a0 - b0;
+  IM(y[2]) = c0 - d0;
+
+  RE(y[3]) = a1 - d1;
+  IM(y[3]) = c1 + b1;  
+}
+
+
+
+
+/*******************************************************************************
+5 points DFT (Winograd algorithm)
+*******************************************************************************/
+void dft5 (complex_t *x,  complex_t* y)
+{
+  complex_t sum[14];
+  complex_t mul[6];
+
+  RE(sum[1]) = RE(x[1]) + RE(x[4]);
+  IM(sum[1]) = IM(x[1]) + IM(x[4]);
+  
+  RE(sum[2]) = RE(x[1]) - RE(x[4]);
+  IM(sum[2]) = IM(x[1]) - IM(x[4]);
+  
+  RE(sum[3]) = RE(x[3]) + RE(x[2]);
+  IM(sum[3]) = IM(x[3]) + IM(x[2]);
+  
+  RE(sum[4]) = RE(x[3]) - RE(x[2]);
+  IM(sum[4]) = IM(x[3]) - IM(x[2]);  
+  
+  RE(sum[5]) = RE(sum[1]) + RE(sum[3]);
+  IM(sum[5]) = IM(sum[1]) + IM(sum[3]);  
+  
+  RE(sum[6]) = RE(sum[1]) - RE(sum[3]);
+  IM(sum[6]) = IM(sum[1]) - IM(sum[3]);  
+  
+  RE(sum[7]) = RE(sum[2]) + RE(sum[4]);
+  IM(sum[7]) = IM(sum[2]) + IM(sum[4]);  
+  
+  RE(y[0]) = RE(sum[5]) + RE(x[0]);
+  IM(y[0]) = IM(sum[5]) + IM(x[0]);
+  
+  RE(mul[1]) = RE(sum[5]) * DFT5_W1;
+  IM(mul[1]) = IM(sum[5]) * DFT5_W1;
+  
+  RE(mul[2]) = RE(sum[6]) * DFT5_W2;
+  IM(mul[2]) = IM(sum[6]) * DFT5_W2;  
+  
+  RE(mul[3]) = -IM(sum[2]) * DFT5_W3;
+  IM(mul[3]) =  RE(sum[2]) * DFT5_W3;
+  
+  RE(mul[4]) = -IM(sum[7]) * DFT5_W4;
+  IM(mul[4]) =  RE(sum[7]) * DFT5_W4;
+  
+  RE(mul[5]) = -IM(sum[4]) * DFT5_W5;
+  IM(mul[5]) =  RE(sum[4]) * DFT5_W5;
+  
+  RE(sum[9]) = RE(y[0]) + RE(mul[1]);
+  IM(sum[9]) = IM(y[0]) + IM(mul[1]);
+  
+  RE(sum[10]) = RE(sum[9]) + RE(mul[2]);
+  IM(sum[10]) = IM(sum[9]) + IM(mul[2]);
+  
+  RE(sum[11]) = RE(sum[9]) - RE(mul[2]);
+  IM(sum[11]) = IM(sum[9]) - IM(mul[2]);
+  
+  RE(sum[12]) = RE(mul[3]) - RE(mul[4]);
+  IM(sum[12]) = IM(mul[3]) - IM(mul[4]);
+  
+  RE(sum[13]) = RE(mul[4]) + RE(mul[5]);
+  IM(sum[13]) = IM(mul[4]) + IM(mul[5]);
+  
+  RE(y[4]) = RE(sum[10]) + RE(sum[12]);
+  IM(y[4]) = IM(sum[10]) + IM(sum[12]);
+  
+  RE(y[3]) = RE(sum[11]) + RE(sum[13]);
+  IM(y[3]) = IM(sum[11]) + IM(sum[13]);
+  
+  RE(y[2]) = RE(sum[11]) - RE(sum[13]);
+  IM(y[2]) = IM(sum[11]) - IM(sum[13]);  
+  
+  RE(y[1]) = RE(sum[10]) - RE(sum[12]);
+  IM(y[1]) = IM(sum[10]) - IM(sum[12]);
+}
+
+
+void dft7 (complex_t *x,  complex_t* y)
+{
+  complex_t sum[31];
+  complex_t mul[9];
+  
+  RE(sum[1]) = RE(x[1]) + RE(x[6]);
+  IM(sum[1]) = IM(x[1]) + IM(x[6]);
+  
+  RE(sum[2]) = RE(x[1]) - RE(x[6]);
+  IM(sum[2]) = IM(x[1]) - IM(x[6]);
+  
+  RE(sum[3]) = RE(x[4]) + RE(x[3]);
+  IM(sum[3]) = IM(x[4]) + IM(x[3]);
+  
+  RE(sum[4]) = RE(x[4]) - RE(x[3]);
+  IM(sum[4]) = IM(x[4]) - IM(x[3]);
+  
+  // Winograd paper mistake?!
+  RE(sum[5]) = RE(x[2]) + RE(x[5]);
+  IM(sum[5]) = IM(x[2]) + IM(x[5]);
+  // Winograd paper mistake?!
+  RE(sum[6]) = RE(x[2]) - RE(x[5]);
+  IM(sum[6]) = IM(x[2]) - IM(x[5]);
+  
+  RE(sum[7]) = RE(sum[1]) + RE(sum[3]);
+  IM(sum[7]) = IM(sum[1]) + IM(sum[3]);
+  
+  RE(sum[8]) = RE(sum[7]) + RE(sum[5]);
+  IM(sum[8]) = IM(sum[7]) + IM(sum[5]);
+  
+  RE(y[0]) = RE(sum[9]) = RE(sum[8]) + RE(x[0]);
+  IM(y[0]) = IM(sum[9]) = IM(sum[8]) + IM(x[0]);
+  
+  RE(sum[10]) = RE(sum[1]) - RE(sum[3]);
+  IM(sum[10]) = IM(sum[1]) - IM(sum[3]);
+  
+  RE(sum[11]) = RE(sum[3]) - RE(sum[5]);
+  IM(sum[11]) = IM(sum[3]) - IM(sum[5]);
+
+  RE(sum[12]) = RE(sum[5]) - RE(sum[1]);
+  IM(sum[12]) = IM(sum[5]) - IM(sum[1]);
+    
+  RE(sum[13]) = RE(sum[2]) + RE(sum[4]);
+  IM(sum[13]) = IM(sum[2]) + IM(sum[4]);
+  
+  RE(sum[14]) = RE(sum[13]) + RE(x[6]);
+  IM(sum[14]) = IM(sum[13]) + IM(x[6]);
+  
+  RE(sum[15]) = RE(sum[2]) - RE(sum[4]);
+  IM(sum[15]) = IM(sum[2]) - IM(sum[4]);
+  
+  RE(sum[16]) = RE(sum[4]) - RE(sum[6]);
+  IM(sum[16]) = IM(sum[4]) - IM(sum[6]);
+  
+  RE(sum[17]) = RE(sum[6]) - RE(sum[2]);
+  IM(sum[17]) = IM(sum[6]) - IM(sum[2]);
+  
+  RE(mul[1]) = DFT7_W1 * RE(sum[8]); 
+  IM(mul[1]) = DFT7_W1 * IM(sum[8]);
+  
+  RE(mul[2]) = DFT7_W2 * RE(sum[10]); 
+  IM(mul[2]) = DFT7_W2 * IM(sum[10]);
+  
+  RE(mul[3]) = DFT7_W3 * RE(sum[11]); 
+  IM(mul[3]) = DFT7_W3 * IM(sum[11]);
+  
+  RE(mul[4]) = DFT7_W4 * RE(sum[12]); 
+  IM(mul[4]) = DFT7_W4 * IM(sum[12]);
+  
+  RE(mul[5]) = -DFT7_W5 * IM(sum[14]); 
+  IM(mul[5]) =  DFT7_W5 * RE(sum[14]);
+  
+  RE(mul[6]) = -DFT7_W6 * IM(sum[15]); 
+  IM(mul[6]) =  DFT7_W6 * RE(sum[15]);
+  
+  RE(mul[7]) = -DFT7_W7 * IM(sum[16]); 
+  IM(mul[7]) =  DFT7_W7 * RE(sum[16]);
+  
+  RE(mul[8]) = -DFT7_W8 * IM(sum[17]); 
+  IM(mul[8]) =  DFT7_W8 * RE(sum[17]);
+  
+  RE(sum[18]) = RE(y[0]) + RE(mul[1]);
+  IM(sum[18]) = IM(y[0]) + IM(mul[1]);
+  
+  RE(sum[19]) = RE(sum[18]) + RE(mul[2]);
+  IM(sum[19]) = IM(sum[18]) + IM(mul[2]);
+  
+  RE(sum[20]) = RE(sum[19]) + RE(mul[3]);
+  IM(sum[20]) = IM(sum[19]) + IM(mul[3]);
+  
+  RE(sum[21]) = RE(sum[18]) - RE(mul[2]);
+  IM(sum[21]) = IM(sum[18]) - IM(mul[2]);
+  
+  RE(sum[22]) = RE(sum[21]) - RE(mul[4]);
+  IM(sum[22]) = IM(sum[21]) - IM(mul[4]);
+  
+  RE(sum[23]) = RE(sum[18]) - RE(mul[3]);
+  IM(sum[23]) = IM(sum[18]) - IM(mul[3]);
+  
+  RE(sum[24]) = RE(sum[23]) + RE(mul[4]);
+  IM(sum[24]) = IM(sum[23]) + IM(mul[4]);
+  
+  RE(sum[25]) = RE(mul[5]) + RE(mul[6]);
+  IM(sum[25]) = IM(mul[5]) + IM(mul[6]);
+  
+  RE(sum[26]) = RE(sum[25]) + RE(mul[7]);
+  IM(sum[26]) = IM(sum[25]) + IM(mul[7]);
+  
+  RE(sum[27]) = RE(mul[5]) - RE(mul[6]);
+  IM(sum[27]) = IM(mul[5]) - IM(mul[6]);
+  
+  RE(sum[28]) = RE(sum[27]) - RE(mul[8]);
+  IM(sum[28]) = IM(sum[27]) - IM(mul[8]);
+  
+  RE(sum[29]) = RE(mul[5]) - RE(mul[7]);
+  IM(sum[29]) = IM(mul[5]) - IM(mul[7]);
+  
+  RE(sum[30]) = RE(sum[29]) + RE(mul[8]);
+  IM(sum[30]) = IM(sum[29]) + IM(mul[8]);
+  
+  RE(y[1]) = RE(sum[20]) + RE(sum[26]);
+  IM(y[1]) = IM(sum[20]) + IM(sum[26]);
+  
+  RE(y[2]) = RE(sum[22]) + RE(sum[28]);
+  IM(y[2]) = IM(sum[22]) + IM(sum[28]);
+  
+  RE(y[3]) = RE(sum[24]) - RE(sum[30]);
+  IM(y[3]) = IM(sum[24]) - IM(sum[30]);
+  
+  RE(y[4]) = RE(sum[24]) + RE(sum[30]);
+  IM(y[4]) = IM(sum[24]) + IM(sum[30]);
+  
+  RE(y[5]) = RE(sum[22]) - RE(sum[28]);
+  IM(y[5]) = IM(sum[22]) - IM(sum[28]);  
+  
+  RE(y[6]) = RE(sum[20]) - RE(sum[26]);
+  IM(y[6]) = IM(sum[20]) - IM(sum[26]);
+  
+}
+
+
+/*******************************************************************************
+16 points DFT (Winograd algorithm)
+*******************************************************************************/
+void dft16 (complex_t *x,  complex_t* y)
+{
+  complex_t t0[16];
+  complex_t t1[16];
+  double tmp;
+  
+  transpose16x16(x, t0);
+  
+  dft4(t0,    t1);
+  dft4(t0+4,  t1+4);
+  dft4(t0+8,  t1+8);
+  dft4(t0+12, t1+12);
+  
+  //#define DFT16_W1       0.923879532511287
+  //#define DFT16_W2       0.382683432365090
+  //#define DFT16_W3       0.707106781186548
+  
+  // 0.923879532511287 - 0.382683432365090i
+  tmp       =  RE(t1[5]) * DFT16_W1 + IM(t1[5]) * DFT16_W2;
+  IM(t1[5]) = -RE(t1[5]) * DFT16_W2 + IM(t1[5]) * DFT16_W1;
+  RE(t1[5]) = tmp;
+  
+  // 0.707106781186548 - 0.707106781186547i
+  tmp       =  ( RE(t1[6]) + IM(t1[6])) * DFT16_W3;
+  IM(t1[6]) =  (-RE(t1[6]) + IM(t1[6])) * DFT16_W3;
+  RE(t1[6]) = tmp;
+    
+  // 0.382683432365090 - 0.923879532511287i
+  tmp       =  RE(t1[7]) * DFT16_W2 + IM(t1[7]) * DFT16_W1;
+  IM(t1[7]) = -RE(t1[7]) * DFT16_W1 + IM(t1[7]) * DFT16_W2;
+  RE(t1[7]) = tmp;
+  
+  // 0.707106781186548 - 0.707106781186547i
+  tmp       =  ( RE(t1[9]) + IM(t1[9])) * DFT16_W3;
+  IM(t1[9]) =  (-RE(t1[9]) + IM(t1[9])) * DFT16_W3;
+  RE(t1[9]) = tmp;
+  
+  // 0.000000000000000 - 1.000000000000000i
+  tmp       =   RE(t1[10]);
+  RE(t1[10]) =  IM(t1[10]);
+  IM(t1[10]) = -tmp;
+  
+  //-0.707106781186547 - 0.707106781186548i
+  tmp        =  (-RE(t1[11]) + IM(t1[11])) * DFT16_W3;
+  IM(t1[11]) =  (-RE(t1[11]) - IM(t1[11])) * DFT16_W3;
+  RE(t1[11]) = tmp;
+  
+  // 0.382683432365090 - 0.923879532511287i
+  tmp       =  RE(t1[13]) * DFT16_W2 + IM(t1[13]) * DFT16_W1;
+  IM(t1[7]) = -RE(t1[13]) * DFT16_W1 + IM(t1[13]) * DFT16_W2;
+  RE(t1[7]) = tmp;
+  
+  //-0.707106781186547 - 0.707106781186548i
+  tmp        =  (-RE(t1[14]) + IM(t1[14])) * DFT16_W3;
+  IM(t1[14]) =  (-RE(t1[14]) - IM(t1[14])) * DFT16_W3;
+  RE(t1[14]) = tmp;
+  
+  
+  //-0.923879532511287 + 0.382683432365090i
+  tmp       =  -RE(t1[15]) * DFT16_W1 - IM(t1[15]) * DFT16_W2;
+  IM(t1[15]) =  RE(t1[15]) * DFT16_W2 - IM(t1[15]) * DFT16_W1;
+  RE(t1[15]) = tmp;
+  
+  transpose16x16(t1, t0);
+  
+  dft4(t0,    t1);
+  dft4(t0+4,  t1+4);
+  dft4(t0+8,  t1+8);
+  dft4(t0+12, t1+12);
+  
+  transpose16x16(t1, y);
+
+}
+
+/*******************************************************************************
+16 x 16 matrix transpose
+*******************************************************************************/
+void transpose16x16(complex_t *x,  complex_t* y)
+{
+  RE(y[ 0]) = RE(x[ 0]);  IM(y[ 0]) = IM(x[ 0]);   
+  RE(y[ 1]) = RE(x[ 4]);  IM(y[ 1]) = IM(x[ 4]);   
+  RE(y[ 2]) = RE(x[ 8]);  IM(y[ 2]) = IM(x[ 8]);   
+  RE(y[ 3]) = RE(x[12]);  IM(y[ 3]) = IM(x[12]);   
+  RE(y[ 4]) = RE(x[ 1]);  IM(y[ 4]) = IM(x[ 1]);   
+  RE(y[ 5]) = RE(x[ 5]);  IM(y[ 5]) = IM(x[ 5]);   
+  RE(y[ 6]) = RE(x[ 9]);  IM(y[ 6]) = IM(x[ 9]);   
+  RE(y[ 7]) = RE(x[13]);  IM(y[ 7]) = IM(x[13]);   
+  RE(y[ 8]) = RE(x[ 2]);  IM(y[ 8]) = IM(x[ 2]);  
+  RE(y[ 9]) = RE(x[ 6]);  IM(y[ 9]) = IM(x[ 6]); 
+  RE(y[10]) = RE(x[10]);  IM(y[10]) = IM(x[10]); 
+  RE(y[11]) = RE(x[14]);  IM(y[11]) = IM(x[14]); 
+  RE(y[12]) = RE(x[ 3]);  IM(y[12]) = IM(x[ 3]); 
+  RE(y[13]) = RE(x[ 7]);  IM(y[13]) = IM(x[ 7]); 
+  RE(y[14]) = RE(x[11]);  IM(y[14]) = IM(x[11]); 
+  RE(y[15]) = RE(x[15]);  IM(y[15]) = IM(x[15]); 
+}
+
+

dspl/src/inout.c

@@ -39,6 +39,78 @@ void DSPL_API dspl_info()
 
 
 
+/******************************************************************************
+Write a real array to the binary file "fn"
+*******************************************************************************/
+int DSPL_API readbin(char* fn, void** x,  int* k, int* dtype)
+{
+  FILE* pFile = NULL;
+  int n, m, t, err;
+  
+  if(!x)
+    return ERROR_PTR;
+  if(!fn)
+    return  ERROR_FNAME;
+
+  pFile = fopen(fn, "rb");
+  if(pFile == NULL)
+      return ERROR_FOPEN;  
+
+  if(fread(&t, sizeof(int), 1, pFile) != 1)
+  {
+    err = ERROR_FREAD_SIZE;
+    goto exit_label;
+  }
+  
+  if(dtype) 
+    *dtype = t;
+  
+  if(fread(&n, sizeof(int), 1, pFile) != 1)
+  {
+    err = ERROR_FREAD_SIZE;
+    goto exit_label;
+  }
+  
+  if(fread(&m, sizeof(int), 1, pFile) != 1)
+  {
+    err = ERROR_FREAD_SIZE;
+    goto exit_label;
+  }
+  
+  if(k)
+    *k = n*m;
+  
+  switch(t)
+  {
+    case DAT_DOUBLE:
+      (*x) = (*x) ? realloc(*x, n*m*sizeof(double)) : 
+                     malloc(n*m*sizeof(double)); 
+      if(fread(*x, sizeof(double), n*m, pFile) != n*m)
+      {
+        err = ERROR_FREAD_SIZE;
+        goto exit_label;
+      }               
+            
+      break;
+    case DAT_COMPLEX:
+      (*x) = (*x) ? realloc(*x, n*m*sizeof(complex_t)) : 
+                     malloc(n*m*sizeof(complex_t)); 
+      if(fread(*x, sizeof(complex_t), n*m, pFile) != n*m)
+      {
+        err = ERROR_FREAD_SIZE;
+        goto exit_label;
+      }   
+    default:
+      err = ERROR_DAT_TYPE;
+      goto exit_label;
+  }
+  err = RES_OK;
+exit_label:
+  if(pFile)
+    fclose(pFile);
+  return err;   
+}
+
 
 
 /******************************************************************************

include/dspl.c

@@ -115,6 +115,7 @@ p_randn                                randn                         ;
 p_randu                                randu                         ;
 p_ratcompos                            ratcompos                     ;
 p_re2cmplx                             re2cmplx                      ;
+p_readbin                              readbin                       ;
 p_signal_pimp                          signal_pimp                   ;
 p_signal_saw                           signal_saw                    ;
 p_sin_cmplx                            sin_cmplx                     ;
@@ -123,6 +124,8 @@ p_sqrt_cmplx                           sqrt_cmplx                    ;
 p_trapint                              trapint                       ;
 p_trapint_cmplx                        trapint_cmplx                 ;
 p_unwrap                               unwrap                        ;
+p_verif                                verif                         ;
+p_verif_cmplx                          verif_cmplx                   ;
 p_writebin                             writebin                      ;
 p_writetxt                             writetxt                      ;
 p_writetxt_3d                          writetxt_3d                   ;
@@ -261,6 +264,7 @@ void* dspl_load()
   LOAD_FUNC(randu);
   LOAD_FUNC(ratcompos);
   LOAD_FUNC(re2cmplx);
+  LOAD_FUNC(readbin);
   LOAD_FUNC(signal_pimp);
   LOAD_FUNC(signal_saw);
   LOAD_FUNC(sin_cmplx);
@@ -269,6 +273,8 @@ void* dspl_load()
   LOAD_FUNC(trapint);
   LOAD_FUNC(trapint_cmplx);
   LOAD_FUNC(unwrap);
+  LOAD_FUNC(verif);
+  LOAD_FUNC(verif_cmplx);
   LOAD_FUNC(writebin);
   LOAD_FUNC(writetxt);
   LOAD_FUNC(writetxt_3d);

include/dspl.h

@@ -100,6 +100,7 @@ typedef struct
 #define ERROR_FILTER_RS                       0x06101819
 #define ERROR_FNAME                           0x06140113
 #define ERROR_FOPEN                           0x06151605
+#define ERROR_FREAD_SIZE                      0x06180501
 #define ERROR_FWRITE_SIZE                     0x06231820
 /* G                                          0x07xxxxxx*/
 /* H                                          0x08xxxxxx*/
@@ -174,7 +175,8 @@ typedef struct
 #define ELLIP_ITER                            16
 #define ELLIP_MAX_ORD                         24
 
-
+#define  DSPL_VERIF_FAILED                    1
+#define  DSPL_VERIF_SUCCESS                   0
 
 #ifdef __cplusplus
   extern "C" {
@@ -648,6 +650,11 @@ DECLARE_FUNC(int,        ratcompos,                   double*          b
 DECLARE_FUNC(int,        re2cmplx,                    double*
                                                 COMMA int
                                                 COMMA complex_t*);
+//------------------------------------------------------------------------------                                                
+DECLARE_FUNC(int,        readbin,                     char*           fn
+                                                COMMA void**          x
+                                                COMMA int*            k
+                                                COMMA int*            dtype);                                                
 //------------------------------------------------------------------------------
 DECLARE_FUNC(int,        signal_pimp,                 double*
                                                 COMMA size_t
@@ -691,6 +698,18 @@ DECLARE_FUNC(int,        unwrap,                      double*
                                                 COMMA int
                                                 COMMA double
                                                 COMMA double);
+//------------------------------------------------------------------------------                                               
+DECLARE_FUNC(int,        verif,                       double*          x
+                                                COMMA double*          y
+                                                COMMA size_t           n
+                                                COMMA double           eps
+                                                COMMA double*          err);    
+//------------------------------------------------------------------------------                                               
+DECLARE_FUNC(int,        verif_cmplx,                 complex_t*       x
+                                                COMMA complex_t*       y
+                                                COMMA size_t           n
+                                                COMMA double           eps
+                                                COMMA double*          err);                                             
 //------------------------------------------------------------------------------
 DECLARE_FUNC(int,        writebin,                    void*
                                                 COMMA int

release/include/dspl.c

@@ -115,6 +115,7 @@ p_randn                                randn                         ;
 p_randu                                randu                         ;
 p_ratcompos                            ratcompos                     ;
 p_re2cmplx                             re2cmplx                      ;
+p_readbin                              readbin                       ;
 p_signal_pimp                          signal_pimp                   ;
 p_signal_saw                           signal_saw                    ;
 p_sin_cmplx                            sin_cmplx                     ;
@@ -123,6 +124,8 @@ p_sqrt_cmplx                           sqrt_cmplx                    ;
 p_trapint                              trapint                       ;
 p_trapint_cmplx                        trapint_cmplx                 ;
 p_unwrap                               unwrap                        ;
+p_verif                                verif                         ;
+p_verif_cmplx                          verif_cmplx                   ;
 p_writebin                             writebin                      ;
 p_writetxt                             writetxt                      ;
 p_writetxt_3d                          writetxt_3d                   ;
@@ -261,6 +264,7 @@ void* dspl_load()
   LOAD_FUNC(randu);
   LOAD_FUNC(ratcompos);
   LOAD_FUNC(re2cmplx);
+  LOAD_FUNC(readbin);
   LOAD_FUNC(signal_pimp);
   LOAD_FUNC(signal_saw);
   LOAD_FUNC(sin_cmplx);
@@ -269,6 +273,8 @@ void* dspl_load()
   LOAD_FUNC(trapint);
   LOAD_FUNC(trapint_cmplx);
   LOAD_FUNC(unwrap);
+  LOAD_FUNC(verif);
+  LOAD_FUNC(verif_cmplx);
   LOAD_FUNC(writebin);
   LOAD_FUNC(writetxt);
   LOAD_FUNC(writetxt_3d);

release/include/dspl.h

@@ -100,6 +100,7 @@ typedef struct
 #define ERROR_FILTER_RS                       0x06101819
 #define ERROR_FNAME                           0x06140113
 #define ERROR_FOPEN                           0x06151605
+#define ERROR_FREAD_SIZE                      0x06180501
 #define ERROR_FWRITE_SIZE                     0x06231820
 /* G                                          0x07xxxxxx*/
 /* H                                          0x08xxxxxx*/
@@ -174,7 +175,8 @@ typedef struct
 #define ELLIP_ITER                            16
 #define ELLIP_MAX_ORD                         24
 
-
+#define  DSPL_VERIF_FAILED                    1
+#define  DSPL_VERIF_SUCCESS                   0
 
 #ifdef __cplusplus
   extern "C" {
@@ -648,6 +650,11 @@ DECLARE_FUNC(int,        ratcompos,                   double*          b
 DECLARE_FUNC(int,        re2cmplx,                    double*
                                                 COMMA int
                                                 COMMA complex_t*);
+//------------------------------------------------------------------------------                                                
+DECLARE_FUNC(int,        readbin,                     char*           fn
+                                                COMMA void**          x
+                                                COMMA int*            k
+                                                COMMA int*            dtype);                                                
 //------------------------------------------------------------------------------
 DECLARE_FUNC(int,        signal_pimp,                 double*
                                                 COMMA size_t
@@ -691,6 +698,18 @@ DECLARE_FUNC(int,        unwrap,                      double*
                                                 COMMA int
                                                 COMMA double
                                                 COMMA double);
+//------------------------------------------------------------------------------                                               
+DECLARE_FUNC(int,        verif,                       double*          x
+                                                COMMA double*          y
+                                                COMMA size_t           n
+                                                COMMA double           eps
+                                                COMMA double*          err);    
+//------------------------------------------------------------------------------                                               
+DECLARE_FUNC(int,        verif_cmplx,                 complex_t*       x
+                                                COMMA complex_t*       y
+                                                COMMA size_t           n
+                                                COMMA double           eps
+                                                COMMA double*          err);                                             
 //------------------------------------------------------------------------------
 DECLARE_FUNC(int,        writebin,                    void*
                                                 COMMA int

test/octave/fft_n1n2.m

@@ -2,9 +2,9 @@ clear all;
 close all;
 clc;
 
-N  = 32;
-N1 = 2;
-N2 = 16;
+N  = 16;
+N1 = 4;
+N2 = 4;
 
 % входной сигнал это вектор столбец размерности [N x 1]
 x = (0:N-1)';

test/octave/fft_twiddle.m

@@ -0,0 +1,22 @@
+clear all; close all; clc;
+
+
+
+
+u = - 2*pi/7;
+
+w = [  (  cos(u) +  cos(2*u) +   cos(3*u)) / 3 - 1;
+       (2*cos(u) -  cos(2*u) -   cos(3*u)) / 3;
+       (  cos(u) -2*cos(2*u) +   cos(3*u)) / 3;
+       (  cos(u) +  cos(2*u) - 2*cos(3*u)) / 3;
+       (  sin(u) +  sin(2*u) -   sin(3*u)) / 3;
+       (2*sin(u) -  sin(2*u) +   sin(3*u)) / 3;
+       (  sin(u) -2*sin(2*u) -   sin(3*u)) / 3;
+       (  sin(u) +  sin(2*u) + 2*sin(3*u)) / 3;]
+   
+   
+fprintf('DFT 7 coeff\n');       
+for k = 1:length(w)
+  fprintf('#define DFT7_W%d      % -.24f\n', k, w(k));
+end  
+         

test/src/fftn_create_test.c

@@ -3,7 +3,7 @@
 #include <string.h>
 #include "dspl.h"
 
-#define N 32
+#define N 7
 int main()
 {
   void* handle;           // DSPL handle
@@ -23,16 +23,9 @@ int main()
   for(n = 0; n < N; n++)
   {
     RE(x[n]) = (double)n;
-    IM(x[n]) = 0.0;
+    IM(x[n]) =0.0;
   }
 
-
-  for(n = 0; n < N*2; n++)
-  {
-    printf("W[%3d] = %12.4f%12.4f\n", n, RE(pfft.w[n]), IM(pfft.w[n]));
-  }
-
-
   fftn_krn(x, y, &pfft, N, 0);
 
 
@@ -41,7 +34,6 @@ int main()
     printf("y[%3d] = %12.4f%12.4f\n", n, RE(y[n]), IM(y[n]));
   }
 
-
   fft_free(&pfft);        // clear FFT object
 
   dspl_free(handle);      // free dspl handle

test/src/fourier_transform_prop_symmetry.c

@@ -0,0 +1,55 @@
+#include <stdio.h>
+#include <string.h>
+#include <math.h>
+#include "dspl.h"
+
+#define N  1000
+
+int main(int argc, char* argv[])
+{
+	double t[N],  s[N], S[N], PHI[N], w[N];
+	complex_t SC[N], sc[N] ;
+	void* handle;
+  int k;
+
+	handle = dspl_load();
+	if(!handle)
+	{
+		printf("cannot to load libdspl!\n");
+		return 0;
+	}
+
+
+	linspace(-1, 10, N, DSPL_SYMMETRIC, t);
+	linspace(-2.0*M_2PI, 2.0*M_2PI, N, DSPL_SYMMETRIC, w);
+
+  for(k =0; k < N; k++)
+  {
+    IM(sc[k]) = 0.0;
+    if(t[k]<0)
+    {
+      s[k] = RE(sc[k]) = 0.0;
+    }
+    else
+    {
+      s[k] = RE(sc[k]) = exp(-t[k]);
+    }
+  }
+
+	writetxt(t, s, N/2, "dat/fourier_transform_prop_sym_time.txt");
+
+  fourier_integral_cmplx(t, sc, N, N, w, SC);
+
+  for(k = 0; k < N; k++)
+  {
+    S[k] = ABS(SC[k]);
+    PHI[k] = ARG(SC[k]);
+  }
+  writetxt(w, S,   N, "dat/fourier_transform_prop_sym_mag.txt");
+  writetxt(w, PHI, N, "dat/fourier_transform_prop_sym_phi.txt");
+
+	// remember to free the resource
+	dspl_free(handle);
+
+	return 0;
+}

verif/.gitignore

@@ -0,0 +1,8 @@
+*.exe
+*.so
+*.dll
+*.exp
+*.ipdb
+*.pdb
+*.iobj
+*.lib

verif/bin/.gitignore

@@ -0,0 +1,8 @@
+*.exe
+*.so
+*.dll
+*.exp
+*.ipdb
+*.pdb
+*.iobj
+*.lib

verif/bin/dat/.gitignore

@@ -0,0 +1,10 @@
+*.dat
+*.txt
+*.exe
+*.so
+*.dll
+*.exp
+*.ipdb
+*.pdb
+*.iobj
+*.lib

verif/obj/.gitignore

@@ -0,0 +1,8 @@
+*.exe
+*.so
+*.dll
+*.exp
+*.ipdb
+*.pdb
+*.iobj
+*.lib

verif/src/readbin_verif.c

@@ -0,0 +1,44 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include "dspl.h"
+
+#define N 10
+int main()
+{
+  void* handle;           // DSPL handle
+  handle = dspl_load();   // Load DSPL function
+  
+  
+  double    *xr = NULL, *yr = NULL;
+  complex_t *xc = NULL, *yc = NULL;
+  double err;
+  int vr;
+  xr = (double*)    malloc(N * sizeof(double));
+  xc = (complex_t*) malloc(N * sizeof(complex_t));
+  
+  randn(xr,   N, 0, 1.0);
+  randn((double*)xc, 2*N, 0, 1.0);
+  
+  writebin(xr, N, DAT_DOUBLE,  "dat/in_real.dat");
+  writebin(xc, N, DAT_COMPLEX, "dat/in_cmplx.dat");
+  
+  readbin("dat/in_real.dat",  (void**)&yr, NULL, NULL);
+  readbin("dat/in_cmplx.dat", (void**)&yc, NULL, NULL);  
+  
+  vr = verif(xr, yr, N, 1E-12, &err);
+  printf("readbin real verification error:          %12.4e\n", err);
+  
+  vr = verif_cmplx(xc, yc, N, 1E-12, &err);
+  printf("readbin cmplx verification error:         %12.4e\n", err);
+  
+  dspl_free(handle);      // free dspl handle
+
+  if(xr) free(xr);
+  if(xc) free(xc);
+  if(yr) free(yr);
+  if(yc) free(yc);
+  return 0;
+}
+
+