Merge pull request #54 from FastNFT/release-0.4

Release 0.4.1

Author: wahls

CHANGELOG.md

@@ -1,5 +1,15 @@
 # Changelog
 
+## [0.4.1] -- 2020-07-13
+
+### Changed
+
+- Number of samples for fnft_nsep again has to be a power of two.
+
+### Fixed
+
+- misc_downsample could return incorrect values for first_last_index[1]
+
 ## [0.4.0] -- 2020-07-08
 
 ### Added

CMakeLists.txt

@@ -21,7 +21,7 @@ project(fnft C)
 
 set(FNFT_VERSION_MAJOR 0)
 set(FNFT_VERSION_MINOR 4)
-set(FNFT_VERSION_PATCH 0)
+set(FNFT_VERSION_PATCH 1)
 set(FNFT_VERSION_SUFFIX "") # should not be longer than FNFT_SUFFIX_MAXLEN
 set(FNFT_VERSION ${FNFT_VERSION_MAJOR}.${FNFT_VERSION_MINOR}.${FNFT_VERSION_PATCH}${FNFT_VERSION_SUFFIX})
 

examples/mex_fnft_nsep_example.m

@@ -21,14 +21,13 @@
 % conditions can be computed using mex_fnft_nsep. The signal is the same
 % plane wave as in the paper https://doi.org/10.1109/TIT.2015.2485944
 
-
 clear all;
 close all;
 
 %%% Setup parameters %%%
 
 T = [0, 2*pi];  % T(1) is the beginning of the period, T(2) is the end
-D = 300;        % number of samples(NOTE: This has to be a even positive integer)
+D = 2^9;        % number of samples (has to be a power of two)
 kappa = +1;     % focusing nonlinear Schroedinger equation  
 
 %%% Setup the signal %%%

include/fnft_nsep.h

@@ -189,7 +189,7 @@ fnft_nsep_opts_t fnft_nsep_default_opts();
  *      - Prins and Wahls, <a href="https://doi.org/10.1109/ICASSP.2018.8461708">&quot;Higher order exponential splittings for the fast non-linear Fourier transform of the KdV equation,&quot;</a>in Proc.ICASSP 2018, Calgary, AB, 2018, pp. 4524-4528.
  *      - Mertsching, <a href="https://doi.org/10.1002/prop.2190350704">&quot; Quasiperiodie Solutions of the Nonlinear Schrödinger Equation,&quot;</a> Fortschr. Phys. 35:519-536, 1987.
  *
- * * The routine supports the following discretizations of type \link fnft_nse_discretization_t \endlink:
+ * The routine supports the following discretizations of type \link fnft_nse_discretization_t \endlink:
  *       - fnft_nse_discretization_2SPLIT1A
  *       - fnft_nse_discretization_2SPLIT1B
  *       - fnft_nse_discretization_2SPLIT2A
@@ -211,7 +211,7 @@ fnft_nsep_opts_t fnft_nsep_default_opts();
  *       - fnft_nse_discretization_2SPLIT8B
  *       - fnft_nse_discretization_4SPLIT4A
  *
- * @param[in] D Number of samples. Has to be even.
+ * @param[in] D Number of samples. Should be power of two and \f$ D\geq 2\f$.
  * @param[in] q Array of length D, contains samples \f$ q(t_n)=q(x_0, t_n) \f$,
  *  where \f$ t_n = T[0] + n*L/D \f$, where \f$L=T[1]-T[0]\f$ is the period and
  *  \f$n=0,1,\dots,D-1\f$, of the to-be-transformed signal in ascending order

matlab/mex_fnft_nsep.c

@@ -64,8 +64,8 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
     kappa = (int)mxGetScalar(prhs[2]);
 
     /* Check values of first three inputs */
-    if ( D<2 || D%2 != 0 )
-        mexErrMsgTxt("Length of the first input q should be even.");
+    if ( D<2 || (D & (D-1)) != 0 )
+        mexErrMsgTxt("Length of the first input q should be >=2 and a power of two.");
     if ( T[0] >= T[1] )
         mexErrMsgTxt("T(1) >= T(2).");
     if ( kappa != +1 && kappa != -1 )
@@ -174,26 +174,26 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
         } else if ( strcmp(str, "quiet") == 0 ) {
 
             fnft_errwarn_setprintf(NULL);
-            
+
         } else if ( strcmp(str, "discr_modal") == 0 ) {
-            
+
             opts.discretization = fnft_nse_discretization_2SPLIT2_MODAL;
-            
+
         } else if ( strcmp(str, "discr_2split2A") == 0 ) {
-            
+
             opts.discretization = fnft_nse_discretization_2SPLIT2A;
-            
+
         } else if ( strcmp(str, "discr_2split4A") == 0 ) {
-            
+
             opts.discretization = fnft_nse_discretization_2SPLIT4A; upsampling_factor = 4;
-            
+
         } else if ( strcmp(str, "discr_2split4B") == 0 ) {
-            
+
             opts.discretization = fnft_nse_discretization_2SPLIT4B; upsampling_factor = 2;
-            
+
         } else if ( strcmp(str, "discr_4split4B") == 0 ) {
-            
-            opts.discretization = fnft_nse_discretization_4SPLIT4B; upsampling_factor = 4;            
+
+            opts.discretization = fnft_nse_discretization_4SPLIT4B; upsampling_factor = 4;
 
         } else {
             snprintf(msg, sizeof msg, "%uth input has invalid value.",

matlab/mex_fnft_nsep.m

@@ -8,8 +8,8 @@
 %   Provides an interface to the C routine fnft_nsep.
 %
 % INPUTS
-%   q               Complex row vector of length D=2n>2. The first sample
-%                   q(1) should correspond to T(1) and the last sample
+%   q               Complex row vector of length D>=2; D should be a power of two.
+%                   The first sample q(1) should correspond to T(1) and the last sample
 %                   q(D) should correspond to T(2)-(T(2)-T(1))/D.
 %   T               Real 1x2 vector with T(2)>T(1).
 %   kappa           +1.0 or -1.0

src/fnft_nsep.c

@@ -92,7 +92,11 @@ INT fnft_nsep(const UINT D, COMPLEX const * const q,
     COMPLEX *q_preprocessed = NULL;
 
     // Check inputs
-    if (D < 2 || D%2 == 1)
+    // D (and Dsub in line 493) only need to be even to ensure that deg in line 279 and 537 are also even for all discretizations. 
+    // deg is required to be even so that its use in lines 321, 355 and 587 make sense. 
+    // Due to the periodic boundary conditions, Dsub has to be an even factor of D. For some choices like D=514, this means Dsub can only be 514 as D=514=2*257.
+    // This makes the idea of subsampling unintuitive. Hence at this stage D is restricted to powers of two for which subsampling still makes sense. 
+    if (D < 2 || (D & (D-1)) != 0 )
         return E_INVALID_ARGUMENT(D);
     if (q == NULL)
         return E_INVALID_ARGUMENT(q);
@@ -483,19 +487,23 @@ static inline INT subsample_and_refine(const UINT D,
             first_last_index, opts_ptr->discretization);
     CHECK_RETCODE(ret_code, release_mem);
 
-    
     // Create a subsampled/resampled version of q for computing initial guesses.
     Dsub = opts_ptr->Dsub;
     if (Dsub == 0) // users wants Dsub to be chosen automatically
-        Dsub = ROUND(SQRT(D * LOG2(D) * LOG2(D)));
-    nskip_per_step = ROUND((REAL)D / Dsub);
-    Dsub = ROUND((REAL)D / nskip_per_step); // actual Dsub
+        Dsub = POW(2.0, CEIL( 0.5 * LOG2(D * LOG2(D) * LOG2(D)) ));
+    else
+        Dsub = POW(2.0, ROUND(LOG2(Dsub))); 
 
     ret_code = nse_discretization_preprocess_signal(D, q, eps_t, kappa, &Dsub, &qsub_preprocessed, &rsub_preprocessed,
             first_last_index, opts_ptr->discretization);
     CHECK_RETCODE(ret_code, release_mem);
+
+    nskip_per_step = D/Dsub;
+    if ( first_last_index[0] != 0 || first_last_index[1]+nskip_per_step != D )
+        return E_ASSERTION_FAILED; // Correct update of T for general
+                                   // downsampling still needs to be
+                                 // implemented
 
-      
     if (upsampling_factor == 2) {        
         nse_discretization = nse_discretization_CF4_2;
     } else if (upsampling_factor == 1) {

test/fnft_nsep/fnft_nsep_test_constant_defocusing_2split4B.c

@@ -25,7 +25,7 @@ INT main()
 {
     INT ret_code;
     const fnft__nsep_testcases_t tc = nsep_testcases_CONSTANT_DEFOCUSING;
-    UINT D = 1026;
+    UINT D = 1024;
     REAL error_bounds[3] = {
         7.6e-5, // main spectrum
         2.5e-8, // aux spectrum

test/fnft_nsep/fnft_nsep_test_constant_defocusing_4split4A.c

@@ -26,10 +26,10 @@ INT main()
 {
     INT ret_code;
     const fnft__nsep_testcases_t tc = nsep_testcases_CONSTANT_DEFOCUSING;
-    UINT D = 68;
+    UINT D = 64;
     REAL error_bounds[3] = {
         6.8e-4, // main spectrum
-        1.2e-6, // aux spectrum
+        1.3e-6, // aux spectrum
         0.0     // sheet indices (zero since not yet implemented)
     };
     fnft_nsep_opts_t opts;