diff --git a/C/README.txt b/C/README.txt new file mode 100644 index 0000000..c9157d9 --- /dev/null +++ b/C/README.txt @@ -0,0 +1,134 @@ + + +C O N T E N T S O F T H E C U R R E N T D I R E C T O R Y + +[./library] - folder with the C code of the AP demodulation library: + + (i) f_ap_demodulation.c - defines the f_ap_demodulation function (the user's + interface to the AP demodulation algorithms). + + (ii) h_ap_demodulation.h - includes all needed external libraries; declares the + input parameter structure for the f_ap_demodulation + function, declares the prototype of the + f_ap_demodulation function itself, defines macros with + error messages and operating-system-dependent new-line + control characters. + + (iii) l_ap_input_validation.c - defines the f_ap_input_validation function, which + is used to check the validity of input parameters + for f_ap_demodulation. + + (iv) l_ap_algorithms.c - defines three functions, namely f_ap_basic, + f_ap_accelerated, and f_ap_projected, which implement + respective AP algorithms. + + (v) l_ap_auxiliary.c - defines the following auxiliary functions used directly + and indirectly by f_ap_demodulation: f_ap_print_error, + f_ap_minmax, f_ap_abs_scaled_max_abs, f_ap_compression, + f_ap_interpolation, f_ap_s_Ub_init, f_ap_mkl_dft_init, + f_ap_mkl_dft_PMw. + +[./examples] - folder with five examples (example[1-5].c) of signal demodulation, + demonstrating various usage cases of f_ap_demodulation. + + + +U S A G E + +The AP demodulation approach was formulated and developed in + + M. Gabrielaitis. "Fast and Accurate Amplitude Demodulation of Wideband Signals," + IEEE Transactions on Signal Processing, vol. 69, pp. 4039-4054, 2021. DOI: + 10.1109/TSP.2021.3087899. + +We refer the user to this publication for any questions related to the working +principle of the new method. + +All demodulation capabilities of the current library are accessible solely via the +f_ap_demodulation function. The sufficient information for the correct use of this +function can be found in the comments section of f_ap_demodulation.c. Specifically, +the requirements for and properties of the input and output arguments are described +there. To consolidate the confidence in using f_ap_demodulation, consider analyzing +the five examples of its applications provided in ./examples. + +To use f_ap_demodulation in an external C code, simply include f_ap_demodulation.c +in the headers of your code files referring to f_ap_demodulation via the +preprocessor's #include directive (see example[1-5].c for an example). The #include +directive with f_ap_demodulation.c can appear in different source files of the same +program if needed - the multiple inclusion error is prevented by an #ifndef...#endif +type guard. For simpler compilation (explained below), consider adding the directory +with the source files of the AP demodulation library to the include path of the +compiler using suitable environment variables. + + + +E X T E R N A L L I B R A R I E S + +Besides the standard C library, f_ap_demodulation uses Intel's oneAPI Math Kernel +Library (oneMKL), specifically, its Fast Fourier Transform routines. Thus, the latter +must be installed on your system. oneMKL is available for all major OS types free of +charge (see https://software.intel.com/content/www/us/en/develop/tools/oneapi/ +base-toolkit/download.html for the installation files). We highly recommend setting +the environment variables with the include and link paths of this library for simpler +compilation - see the installation guide of the MKL library for the instructions. + + + +C O M P I L A T I O N + +C code using the AP demodulation library is supposed to be compiled with GCC. This is +the state-of-the-art compiler available free of charge for all major OS types. It is +preinstalled on Linux systems as a rule. A GCC installation guide for Windows and Mac +can be found following this link: https://www.guru99.com/c-gcc-install.html. + +The exact directive to the compiler depends on whether the OS environment variables +corresponding to the include, link, and load paths of the compiler are set or not: + +(i) If no include path of the AP demodulation library and no include, link, and load + paths of the oneMKL library are set via appropriate environment variables, the + basic form of the compilation directive used to generate an executable file is + + gcc -I Path_to_AP_Demodulation -I Path_to_MKL_include mycode.c \ + Path_to_MKL_lib/libmkl_rt.so -lm -o myprogram + + Here, + + "Path_to_AP_Demodulation" is the full path of the folder with the C code of + the AP demodulation library, e.g., /usr/local/ap-demodulation/C/library/ + + "Path_to_MKL_include" is the full path of the folder with the header files of + the MKL library. For example, if a 2019 release of the MKL was installed in + /usr/local/intel on a Linux system, Path_to_MKL_include would be + /usr/local/intel/compilers_and_libraries_2019/linux/mkl/include + + "mycode.c" is the C code that uses f_ap_demodulation to be compiled. + + "Path_to_MKL_lib/libmkl_rt.so" is the full path of the needed MKL (shared) + library file that contains precompiled FFT routines. For example, if a 2019 + release of the MKL was installed in /usr/local/intel on a Linux system, + Path_to_MKL_lib would be + /usr/local/intel/compilers_and_libraries_2019/linux/mkl/lib/intel64 + + "-lm" refers to the full path to the standard C (shared) library file that + contains precompiled basic math routines. + + "myprogram" is the full path of the executable to be generated. + +(ii) If the include and load paths of the AP demodulation and oneMKL libraries are + set via the environment variables, the compiler directive for generating an + executable file simplifies to + + gcc mycode.c -lmkl_rt -lm -o myprogram + +Note that, in both cases, the MKL library is dynamically linked. That is, myprogram +cannot run successfully without relevant MKL library files present on the system +during runtime. + +The above compilation statements are minimal in the sense that they do not assume +additional source files, libraries, or compilation options. The latter can be added +following the standard input argument syntax of GCC if needed. Please check +"An Introduction to GCC" by Brian Gough, an excellent and concise text on the most +relevant aspects of GCC usage in practice. In particular, Section 3.1 there explains +how to set the environment variables relevant for the compilation process on a Linux +system. In the case of another OS, google for analogs. + diff --git a/C/examples/example1.c b/C/examples/example1.c new file mode 100644 index 0000000..3d9bd57 --- /dev/null +++ b/C/examples/example1.c @@ -0,0 +1,349 @@ + +/* C O P Y R I G H T N O T I C E + * + * Copyright ©2021. Institute of Science and Technology Austria (IST Austria). + * All Rights Reserved. The underlying technology is protected by PCT Patent + * Application No. PCT/EP2021/054650. + * + * This file is part of the AP demodulation library, which is free software: you can + * redistribute it and/or modify it under the terms of the GNU General Public License + * as published by the Free Software Foundation in version 2. + * + * This program 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 v2 for more details. You + * should have received a copy of the GNU General Public License v2 along with this + * program. If not, see https://www.gnu.org/licenses/. + * + * Contact the Technology Transfer Office, IST Austria, Am Campus 1, + * A-3400 Klosterneuburg, Austria, +43-(0)2243 9000, twist@ist.ac.at, for commercial + * licensing opportunities. + * + * All other inquiries should be directed to the author, Mantas Gabrielaitis, + * mgabriel@ist.ac.at + */ + + +/* EXAMPLE 1 + * + * In this example, a synthetic 1D amplitude-modulated signal built from a harmonic + * carrier and a sinusoidal modulator is generated and demodulated by using the + * AP-Basic algorithm. Sample points of the the predefined and inferred modulators + * and carriers are then written into a text file for further analysis. This example + * illustrates application of the function 'f_ap_demodulation' in its simplest form, + * when the signal is sampled uniformly, and when only the final estimates of the + * modulator and infeasibility error are required. */ + + + +#include + +#include + +#include "f_ap_demodulation.c" + + + +#ifdef _WIN32 + + #define STR_NL "\r\n" + +#else + + #define STR_NL "\n" + +#endif + + + +int main(void) +{ + + /* Declarations and initializations */ + + int line; + + int exitflag = 0; + + + long i, j, iter, n; + + double dt; + + double *m=NULL, *c=NULL, *s=NULL, *Ub=NULL, *t=NULL; + + double *out_m=NULL, *out_c=NULL, *out_e=NULL; + + double w[] = {0.4170, 0.7203, 0.0001, 0.3023, 0.1468, 0.0923, 0.1863, 0.3456, \ + 0.3968, 0.5388, 0.4192, 0.6852, 0.2045, 0.8781, 0.0274, 0.6705, \ + 0.4173, 0.5587, 0.1404, 0.1981, 0.8007, 0.9683, 0.3134, 0.6923, \ + 0.8764, 0.8946, 0.0850, 0.0391, 0.1698, 0.8781, 0.0983, 0.4211, \ + 0.9579, 0.5332, 0.6919, 0.3155, 0.6865, 0.8346, 0.0183, 0.7501}; + + struct s_ParamsAP Par; + + Par.Fs = NULL; Par.Fc = NULL; Par.Ns = NULL; Par.im = NULL; Par.ie = NULL; + + FILE *fid = NULL; + + + + /* Number of sample points */ + + n = 16384; + + + + /* Time step */ + + dt = (double)25/(n-1); + + + + /* Modulator (nonnegative sinusoidal) */ + + m = (double*) malloc(n*sizeof(double)); + + if (m == NULL) {line=__LINE__-2; exitflag=-1; goto failed;} + + + for (i=0; i + +#include + +#include "f_ap_demodulation.c" + + + +#ifdef _WIN32 + + #define STR_NL "\r\n" + +#else + + #define STR_NL "\n" + +#endif + + + +int main(void) +{ + + /* Declarations and initializations */ + + int line; + + int exitflag = 0; + + + long i0, i1, j, i_lin, iter, n[2]; + + double dt[2]; + + double *m=NULL, *c=NULL, *s=NULL, *Ub=NULL, *t=NULL; + + double *out_m=NULL, *out_c=NULL, *out_e=NULL; + + double w[] = {0.5173, 0.9470, 0.7655, 0.2824, 0.2210, 0.6862, 0.1671, 0.3924, \ + 0.6181, 0.4119, 0.0025, 0.8840}; + + double cc[] = {0.2220, 0.2067, 0.4884, 0.7659, 0.2968, 0.0807, 0.4413, 0.8799, \ + 0.4142, 0.6288, 0.5999, 0.2847, 0.3276, 0.1656, 0.9602, 0.0243, \ + 0.6998, 0.0229, 0.0016, 0.6398, 0.2591, 0.8705, 0.0022, 0.9815, \ + 0.8137, 0.0291, 0.1115, 0.9649, 0.6354, 0.9267, 0.8248, 0.3610, \ + 0.5464, 0.3655, 0.7951, 0.6389, 0.5835, 0.9435, 0.8436, 0.1008, \ + 0.5104, 0.4783, 0.5147, 0.8005, 0.5726, 0.9851, 0.4524, 0.3320, \ + 0.4077, 0.3303, 0.5267, 0.8930, 0.7699, 0.7100, 0.7673, 0.2396, \ + 0.3636, 0.0601, 0.8132, 0.0634, 0.0851, 0.1703, 0.8146, 0.0598, \ + 0.5710, 0.8257, 0.5809, 0.6523, 0.4593, 0.4967, 0.6040, 0.2185, \ + 0.3653, 0.5421, 0.1344, 0.6792, 0.9972, 0.7268, 0.5885, 0.2456, \ + 0.1656, 0.2782, 0.6248, 0.8938, 0.7703, 0.1198, 0.1547, 0.0004, \ + 0.1611, 0.3383, 0.6337, 0.7140, 0.1907, 0.9386, 0.2925, 0.5712, \ + 0.5882, 0.6689, 0.2065, 0.8112, 0.4473, 0.0477, 0.7143, 0.9132, \ + 0.5831, 0.2310, 0.5443, 0.0688, 0.2748, 0.4697, 0.9382, 0.7194, \ + 0.6532, 0.4846, 0.2203, 0.7646, 0.3682, 0.2952, 0.3210, 0.8742, \ + 0.3407, 0.3807, 0.8874, 0.0466, 0.4245, 0.6895, 0.8714, 0.1656, \ + 0.2115, 0.3212, 0.0040, 0.5278, 0.2544, 0.0021, 0.7132, 0.2993, \ + 0.1577, 0.0797, 0.4320, 0.7166, 0.6996, 0.2169, 0.4249, 0.3518, \ + 0.9538, 0.0420, 0.2370, 0.2778, 0.1176, 0.9295, 0.9049, 0.9100, \ + 0.0045, 0.4728, 0.2697, 0.8605, 0.3467, 0.1765, 0.1566, 0.3896, \ + 0.3706, 0.0925, 0.2739, 0.5697, 0.1154, 0.5252, 0.1085, 0.6850, \ + 0.5321, 0.2710, 0.4101, 0.4250, 0.6209, 0.4344, 0.0603, 0.6348, \ + 0.8840, 0.1613, 0.2017, 0.9767, 0.4458, 0.0140, 0.5626, 0.0921, \ + 0.3696, 0.1786, 0.5415, 0.3791, 0.4891, 0.9199, 0.7476, 0.0325, \ + 0.8198, 0.6506, 0.4079, 0.5986, 0.1102, 0.9745, 0.2820, 0.5881, \ + 0.6075, 0.0562, 0.1923, 0.9325, 0.5379, 0.8779, 0.8283, 0.9809, \ + 0.7651, 0.7090, 0.9867, 0.8185, 0.5057, 0.9999, 0.2444, 0.0221, \ + 0.4203, 0.9045, 0.7554, 0.7669, 0.2311, 0.5881, 0.3348, 0.6998, \ + 0.6015, 0.8316, 0.9592, 0.6635, 0.7464, 0.5475, 0.4010, 0.8478, \ + 0.7038, 0.6468, 0.2885, 0.3513, 0.8803, 0.4530, 0.8347, 0.8748, \ + 0.1201, 0.4327, 0.1652, 0.9487, 0.1985, 0.1338, 0.5014, 0.0863, \ + 0.9117, 0.2444, 0.0224, 0.7622, 0.1201, 0.5540, 0.4990, 0.0593, \ + 0.4552, 0.8290, 0.7764, 0.3111, 0.0139, 0.4933, 0.9800, 0.9281, \ + 0.5271, 0.6275, 0.8195, 0.2567, 0.6348, 0.9209, 0.8707, 0.9186, \ + 0.6117, 0.5184, 0.1877, 0.7384, 0.1583, 0.2741, 0.2961, 0.5798, \ + 0.2658, 0.2536, 0.1442, 0.9639, 0.1884, 0.2046, 0.7795, 0.5777, \ + 0.5155, 0.9856, 0.8025, 0.9227, 0.4695, 0.3989, 0.5465, 0.0865, \ + 0.2512, 0.6318, 0.8119, 0.9126, 0.0942, 0.0355, 0.7961, 0.2443, \ + 0.3521, 0.4934, 0.9393, 0.1117, 0.3460, 0.3834, 0.9611, 0.8500, \ + 0.4466, 0.0204, 0.4114, 0.8014, 0.7029, 0.3600, 0.8986, 0.0827, \ + 0.6608, 0.9652, 0.7591, 0.7016, 0.9743, 0.2215, 0.8103, 0.4497, \ + 0.2642, 0.2421, 0.8078, 0.1953, 0.8103, 0.9650, 0.3025, 0.6594, \ + 0.4728, 0.2419, 0.4187, 0.2111, 0.7486, 0.4256, 0.4779, 0.2652, \ + 0.3019, 0.6013, 0.5609, 0.6158, 0.6050, 0.0713, 0.0604, 0.3483, \ + 0.1597, 0.7306, 0.1815, 0.9303, 0.4071, 0.8455, 0.7393, 0.7428, \ + 0.8889, 0.2651, 0.6342, 0.2582, 0.4256, 0.9901, 0.8879, 0.3714, \ + 0.3307, 0.4458, 0.0275, 0.8636, 0.8304, 0.4334, 0.3488, 0.0260, \ + 0.0195, 0.5254, 0.9706, 0.2852, 0.9667, 0.1369, 0.7670, 0.0974, \ + 0.7489, 0.6573, 0.6528, 0.7640, 0.5043, 0.4854, 0.2949, 0.5953, \ + 0.0584, 0.3418, 0.8476, 0.4459, 0.1661, 0.5777, 0.9752, 0.0141, \ + 0.7735, 0.7273, 0.9084, 0.6823, 0.1977, 0.8807, 0.8074, 0.0984, \ + 0.5922, 0.5683, 0.3182, 0.4868, 0.0836, 0.6991, 0.0567, 0.6351, \ + 0.7000, 0.6762, 0.1547, 0.6705, 0.6147, 0.3651, 0.2509, 0.9189, \ + 0.4966, 0.5290, 0.8573, 0.8691, 0.2597, 0.7464, 0.4243, 0.0416, \ + 0.1331, 0.8380, 0.6698, 0.3229, 0.8539, 0.9468, 0.9107, 0.5435, \ + 0.5347, 0.9536, 0.4167, 0.4788, 0.1567, 0.6897, 0.5744, 0.5319, \ + 0.2763, 0.7413, 0.2140, 0.8967, 0.1410, 0.5868, 0.1654, 0.1175, \ + 0.0989, 0.1180, 0.5080, 0.8384, 0.3963, 0.0145, 0.9733, 0.5382, \ + 0.3626, 0.5241, 0.4665, 0.8944, 0.0441, 0.0947, 0.2165, 0.3988, \ + 0.6922, 0.5553, 0.1900, 0.6840, 0.6099, 0.2820, 0.8589, 0.5188, \ + 0.6836, 0.1999, 0.7759, 0.1299, 0.6388, 0.3019, 0.6396, 0.0041, \ + 0.0733, 0.2318, 0.5796, 0.0011, 0.9319, 0.8927, 0.6560, 0.7422, \ + 0.3567, 0.0286, 0.0101, 0.5152, 0.6869, 0.0101, 0.0947, 0.3179, \ + 0.6055, 0.9869, 0.8656, 0.3513, 0.8832, 0.8192, 0.9975, 0.7987, \ + 0.4752, 0.4678, 0.3848, 0.1052, 0.5652, 0.3117, 0.5121, 0.1686, \ + 0.7248, 0.9140, 0.7688, 0.5681, 0.3609, 0.6887, 0.9962, 0.4963, \ + 0.6052, 0.5642, 0.1326, 0.1580, 0.0161, 0.7922, 0.4014, 0.0568, \ + 0.3983, 0.3021, 0.1127, 0.3214}; + + struct s_ParamsAP Par; + + Par.Fs = NULL; Par.Fc = NULL; Par.Ns = NULL; Par.im = NULL; Par.ie = NULL; + + FILE *fid; + + + + /* Number of sample points in each dimension */ + + n[0] = 201; + + n[1] = 201; + + + + /* Time steps */ + + dt[0] = 0.005; + + dt[1] = 0.005; + + + + /* Modulator (LP-random-field) */ + + m = (double*) calloc(n[0]*n[1], sizeof(double)); + + if (m == NULL) {line=__LINE__-2; exitflag=-1; goto failed;} + + + for (i0=0; i0 + +#include + +#include "f_ap_demodulation.c" + + + +#ifdef _WIN32 + + #define STR_NL "\r\n" + +#else + + #define STR_NL "\n" + +#endif + + + +int main(void) +{ + + /* Declarations and initializations */ + + int line; + + int exitflag = 0; + + + long i, j, iter, n; + + double *m=NULL, *c=NULL, *s=NULL, *Ub=NULL; + + double *out_m=NULL, *out_c=NULL, *out_e=NULL; + + double t[] = {0.0000, 0.0087, 0.0213, 0.0246, 0.0318, 0.0370, 0.0414, 0.0471, \ + 0.0549, 0.0633, 0.0736, 0.0823, 0.0945, 0.1004, 0.1151, 0.1187, \ + 0.1307, 0.1394, 0.1500, 0.1551, 0.1609, 0.1746, 0.1905, 0.1978, \ + 0.2101, 0.2248, 0.2397, 0.2441, 0.2478, 0.2533, 0.2680, 0.2725, \ + 0.2813, 0.2970, 0.3072, 0.3195, 0.3269, 0.3391, 0.3532, 0.3567, \ + 0.3698, 0.3859, 0.3989, 0.4058, 0.4194, 0.4240, 0.4331, 0.4482, \ + 0.4553, 0.4623, 0.4672, 0.4707, 0.4828, 0.4888, 0.4956, 0.5052, \ + 0.5092, 0.5199, 0.5251, 0.5360, 0.5484, 0.5530, 0.5617, 0.5740, \ + 0.5826, 0.5865, 0.5968, 0.6087, 0.6187, 0.6342, 0.6451, 0.6602, \ + 0.6652, 0.6703, 0.6841, 0.6925, 0.6979, 0.7133, 0.7211, 0.7341, \ + 0.7468, 0.7616, 0.7730, 0.7860, 0.7938, 0.8006, 0.8155, 0.8244, \ + 0.8402, 0.8521, 0.8635, 0.8682, 0.8839, 0.8930, 0.9038, 0.9124, \ + 0.9187, 0.9337, 0.9445, 0.9478, 0.9591, 0.9666, 0.9767, 0.9915, \ + 0.9994, 1.0145, 1.0259, 1.0294, 1.0447, 1.0570, 1.0733, 1.0788, \ + 1.0838, 1.0992, 1.1116, 1.1157, 1.1288, 1.1419, 1.1572, 1.1697, \ + 1.1746, 1.1781, 1.1817, 1.1853, 1.1918, 1.2063, 1.2165, 1.2270, \ + 1.2412, 1.2461, 1.2530, 1.2639, 1.2798, 1.2904, 1.2939, 1.3076, \ + 1.3138, 1.3276, 1.3359, 1.3505, 1.3634, 1.3740, 1.3790, 1.3830, \ + 1.3879, 1.3917, 1.3964, 1.4026, 1.4151, 1.4257, 1.4291, 1.4333, \ + 1.4492, 1.4598, 1.4657, 1.4723, 1.4852, 1.4910, 1.5019, 1.5178, \ + 1.5321, 1.5384, 1.5481, 1.5595, 1.5735, 1.5788, 1.5823, 1.5865, \ + 1.5961, 1.6073, 1.6179, 1.6253, 1.6415, 1.6523, 1.6605, 1.6709, \ + 1.6839, 1.6959, 1.7026, 1.7067, 1.7148, 1.7263, 1.7323, 1.7453, \ + 1.7495, 1.7561, 1.7699, 1.7756, 1.7872, 1.7973, 1.8126, 1.8193, \ + 1.8235, 1.8363, 1.8496, 1.8647, 1.8801, 1.8835, 1.8899, 1.9012, \ + 1.9168, 1.9324, 1.9429, 1.9581, 1.9697, 1.9781, 1.9877, 1.9988, \ + 2.0092, 2.0246, 2.0398, 2.0482, 2.0640, 2.0695, 2.0744, 2.0795, \ + 2.0893, 2.0929, 2.1085, 2.1225, 2.1260, 2.1315, 2.1391, 2.1441, \ + 2.1579, 2.1656, 2.1811, 2.1920, 2.2067, 2.2210, 2.2360, 2.2453, \ + 2.2556, 2.2693, 2.2763, 2.2859, 2.2970, 2.3005, 2.3115, 2.3204, \ + 2.3342, 2.3415, 2.3564, 2.3672, 2.3729, 2.3864, 2.3976, 2.4016, \ + 2.4103, 2.4224, 2.4377, 2.4409, 2.4569, 2.4651, 2.4810, 2.4922, \ + 2.5062, 2.5170, 2.5284, 2.5354, 2.5463, 2.5593, 2.5738, 2.5869, \ + 2.5992, 2.6138, 2.6212, 2.6332, 2.6424, 2.6506, 2.6592, 2.6677, \ + 2.6751, 2.6865, 2.6953, 2.7113, 2.7234, 2.7292, 2.7380, 2.7458, \ + 2.7594, 2.7741, 2.7892, 2.8011, 2.8079, 2.8144, 2.8288, 2.8389, \ + 2.8527, 2.8634, 2.8762, 2.8862, 2.8995, 2.9102, 2.9195, 2.9272, \ + 2.9314, 2.9396, 2.9439, 2.9599, 2.9656, 2.9794, 2.9941, 3.0063, \ + 3.0170, 3.0223, 3.0317, 3.0394, 3.0456, 3.0566, 3.0639, 3.0791, \ + 3.0942, 3.1009, 3.1056, 3.1113, 3.1211, 3.1339, 3.1398, 3.1463, \ + 3.1607, 3.1693, 3.1806, 3.1869, 3.1915, 3.2015, 3.2110, 3.2162, \ + 3.2276, 3.2380, 3.2497, 3.2549, 3.2679, 3.2741, 3.2841, 3.2976, \ + 3.3012, 3.3086, 3.3233, 3.3375, 3.3478, 3.3624, 3.3780, 3.3920, \ + 3.4064, 3.4110, 3.4227, 3.4351, 3.4464, 3.4600, 3.4637, 3.4770, \ + 3.4898, 3.4965, 3.5031, 3.5146, 3.5223, 3.5360, 3.5451, 3.5585, \ + 3.5747, 3.5819, 3.5870, 3.6020, 3.6123, 3.6283, 3.6398, 3.6560, \ + 3.6664, 3.6765, 3.6815, 3.6894, 3.6930, 3.6984, 3.7114, 3.7150, \ + 3.7231, 3.7376, 3.7498, 3.7621, 3.7678, 3.7768, 3.7877, 3.8038, \ + 3.8097, 3.8162, 3.8229, 3.8359, 3.8452, 3.8492, 3.8590, 3.8651, \ + 3.8787, 3.8859, 3.8895, 3.9005, 3.9147, 3.9230, 3.9360, 3.9459, \ + 3.9562, 3.9720, 3.9857, 3.9894, 4.0019, 4.0112, 4.0268, 4.0330, \ + 4.0397, 4.0440, 4.0529, 4.0575, 4.0691, 4.0828, 4.0951, 4.1084, \ + 4.1161, 4.1304, 4.1392, 4.1532, 4.1646, 4.1698, 4.1741, 4.1775, \ + 4.1817, 4.1909, 4.1956, 4.1993, 4.2124, 4.2208, 4.2338, 4.2429, \ + 4.2520, 4.2615, 4.2710, 4.2847, 4.2932, 4.3082, 4.3120, 4.3253, \ + 4.3302, 4.3415, 4.3449, 4.3552, 4.3585, 4.3742, 4.3892, 4.4029, \ + 4.4180, 4.4232, 4.4285, 4.4342, 4.4456, 4.4607, 4.4768, 4.4893, \ + 4.5021, 4.5172, 4.5257, 4.5322, 4.5378, 4.5426, 4.5564, 4.5616, \ + 4.5683, 4.5822, 4.5895, 4.6056, 4.6123, 4.6225, 4.6299, 4.6450, \ + 4.6531, 4.6620, 4.6719, 4.6874, 4.6911, 4.7037, 4.7185, 4.7222, \ + 4.7322, 4.7397, 4.7542, 4.7647, 4.7770, 4.7861, 4.7976, 4.8046, \ + 4.8080, 4.8188, 4.8261, 4.8361, 4.8513, 4.8601, 4.8666, 4.8747, \ + 4.8901, 4.9056, 4.9198, 4.9351, 4.9413, 4.9457, 4.9519, 4.9593, \ + 4.9648, 4.9760, 4.9846, 4.9985, 5.0042, 5.0166, 5.0230, 5.0337, \ + 5.0415, 5.0455, 5.0518, 5.0637, 5.0734, 5.0835, 5.0890, 5.0997, \ + 5.1159, 5.1298, 5.1409, 5.1568, 5.1718, 5.1829, 5.1865, 5.1910, \ + 5.1951, 5.2043, 5.2124, 5.2284, 5.2338, 5.2498, 5.2630, 5.2770, \ + 5.2885, 5.3005, 5.3100, 5.3134, 5.3213, 5.3310, 5.3437, 5.3531, \ + 5.3623, 5.3674, 5.3708, 5.3839, 5.3913, 5.4074, 5.4135, 5.4212, \ + 5.4313, 5.4444, 5.4537, 5.4586, 5.4659, 5.4758, 5.4878, 5.5011, \ + 5.5061, 5.5096, 5.5196, 5.5334, 5.5369, 5.5489, 5.5611, 5.5702, \ + 5.5854, 5.5970, 5.6004, 5.6099, 5.6244, 5.6385, 5.6502, 5.6622, \ + 5.6730, 5.6799, 5.6904, 5.7024, 5.7103, 5.7247, 5.7305, 5.7435, \ + 5.7505, 5.7639, 5.7727, 5.7865, 5.7944, 5.8004, 5.8137, 5.8210, \ + 5.8338, 5.8467, 5.8529, 5.8589, 5.8648, 5.8699, 5.8781, 5.8817, \ + 5.8864, 5.8984, 5.9121, 5.9164, 5.9227, 5.9287, 5.9439, 5.9564, \ + 5.9669, 5.9741, 5.9882, 5.9972, 6.0125, 6.0249, 6.0344, 6.0393, \ + 6.0553, 6.0606, 6.0665, 6.0754, 6.0839, 6.0891, 6.1019, 6.1076, \ + 6.1192, 6.1323, 6.1383, 6.1494, 6.1624, 6.1740, 6.1851, 6.1922, \ + 6.2050, 6.2205, 6.2293, 6.2428, 6.2468, 6.2609, 6.2667, 6.2751, \ + 6.2823, 6.2866, 6.3016, 6.3097, 6.3199, 6.3296, 6.3346, 6.3405, \ + 6.3448, 6.3546, 6.3613, 6.3692, 6.3739, 6.3874, 6.3921, 6.4082, \ + 6.4137, 6.4244, 6.4283, 6.4418, 6.4475, 6.4577, 6.4706, 6.4758, \ + 6.4862, 6.4923, 6.5055, 6.5181, 6.5237, 6.5382, 6.5417, 6.5562, \ + 6.5667, 6.5752, 6.5884, 6.6010, 6.6171, 6.6240, 6.6273, 6.6427, \ + 6.6572, 6.6699, 6.6799, 6.6924, 6.7058, 6.7140, 6.7273, 6.7403, \ + 6.7516, 6.7600, 6.7724, 6.7757, 6.7891, 6.8040, 6.8104, 6.8152, \ + 6.8213, 6.8285, 6.8433, 6.8536, 6.8606, 6.8657, 6.8727, 6.8840, \ + 6.8915, 6.9007, 6.9097, 6.9238, 6.9326, 6.9404, 6.9524, 6.9586, \ + 6.9679, 6.9753, 6.9867, 7.0014, 7.0105, 7.0240, 7.0332, 7.0450, \ + 7.0500, 7.0589, 7.0740, 7.0851, 7.0984, 7.1082, 7.1180, 7.1322, \ + 7.1364, 7.1471, 7.1626, 7.1726, 7.1784, 7.1927, 7.1993, 7.2117, \ + 7.2220, 7.2376, 7.2490, 7.2632, 7.2665, 7.2827, 7.2869, 7.2944, \ + 7.3100, 7.3134, 7.3274, 7.3418, 7.3508, 7.3574, 7.3711, 7.3806, \ + 7.3856, 7.4010, 7.4159, 7.4256, 7.4400, 7.4487, 7.4609, 7.4693, \ + 7.4792, 7.4849, 7.5007, 7.5078, 7.5124, 7.5176, 7.5210, 7.5336, \ + 7.5442, 7.5578, 7.5677, 7.5813, 7.5936, 7.6070, 7.6156, 7.6273, \ + 7.6329, 7.6403, 7.6458, 7.6544, 7.6608, 7.6694, 7.6853, 7.6928, \ + 7.7088, 7.7204, 7.7285, 7.7430, 7.7543, 7.7608, 7.7744, 7.7833, \ + 7.7913, 7.7988, 7.8112, 7.8179, 7.8317, 7.8388, 7.8492, 7.8588, \ + 7.8732, 7.8880, 7.8937, 7.9046, 7.9196, 7.9286, 7.9439, 7.9508, \ + 7.9620, 7.9741, 7.9804, 7.9838, 7.9925, 8.0080, 8.0157, 8.0291, \ + 8.0347, 8.0424, 8.0475, 8.0601, 8.0725, 8.0847, 8.0913, 8.1036, \ + 8.1098, 8.1186, 8.1267, 8.1346, 8.1386, 8.1501, 8.1626, 8.1738, \ + 8.1855, 8.1910, 8.1962, 8.2062, 8.2208, 8.2265, 8.2358, 8.2446, \ + 8.2544, 8.2597, 8.2674, 8.2741, 8.2884, 8.3021, 8.3109, 8.3221, \ + 8.3272, 8.3371, 8.3442, 8.3587, 8.3707, 8.3822, 8.3871, 8.3965, \ + 8.4126, 8.4282, 8.4399, 8.4451, 8.4567, 8.4674, 8.4767, 8.4856, \ + 8.4966, 8.5110, 8.5240, 8.5348, 8.5501, 8.5542, 8.5704, 8.5744, \ + 8.5802, 8.5890, 8.5936, 8.6050, 8.6089, 8.6159, 8.6199, 8.6324, \ + 8.6443, 8.6525, 8.6582, 8.6712, 8.6789, 8.6925, 8.7022, 8.7123, \ + 8.7159, 8.7276, 8.7354, 8.7416, 8.7505, 8.7588, 8.7682, 8.7842, \ + 8.7922, 8.8056, 8.8160, 8.8309, 8.8387, 8.8452, 8.8603, 8.8642, \ + 8.8798, 8.8903, 8.8985, 8.9147, 8.9187, 8.9287, 8.9324, 8.9431, \ + 8.9487, 8.9602, 8.9762, 8.9909, 9.0000, 9.0125, 9.0259, 9.0356, \ + 9.0458, 9.0510, 9.0591, 9.0642, 9.0769, 9.0864, 9.0955, 9.1103, \ + 9.1204, 9.1290, 9.1358, 9.1400, 9.1487, 9.1523, 9.1593, 9.1692, \ + 9.1850, 9.1897, 9.2017, 9.2115, 9.2249, 9.2300, 9.2343, 9.2428, \ + 9.2565, 9.2622, 9.2755, 9.2825, 9.2886, 9.2921, 9.3005, 9.3088, \ + 9.3206, 9.3248, 9.3300, 9.3335, 9.3383, 9.3500, 9.3585, 9.3660, \ + 9.3765, 9.3927, 9.4068, 9.4192, 9.4344, 9.4382, 9.4424, 9.4518, \ + 9.4596, 9.4751, 9.4847, 9.4950, 9.5100, 9.5190, 9.5337, 9.5403, \ + 9.5471, 9.5546, 9.5650, 9.5712, 9.5832, 9.5883, 9.5928, 9.6074, \ + 9.6137, 9.6220, 9.6327, 9.6428, 9.6471, 9.6617, 9.6774, 9.6912, \ + 9.6982, 9.7083, 9.7160, 9.7265, 9.7425, 9.7498, 9.7618, 9.7693, \ + 9.7826, 9.7901, 9.8050, 9.8180, 9.8312, 9.8406, 9.8466, 9.8504, \ + 9.8579, 9.8616, 9.8739, 9.8859, 9.8948, 9.9080, 9.9183, 9.9221, \ + 9.9271, 9.9328, 9.9405, 9.9444, 9.9556, 9.9655, 9.9768, 9.9857}; + + double w[] = {1.5648, 0.5312, 0.1413, 0.7588, -0.8616, -0.3586, 0.9106, -0.1787,\ + -0.0108, -0.0989, -0.3559, -0.4015, 0.2917, -0.3458, -1.1990, \ + 0.7651, -0.9884, -1.1668, 0.6584, -1.3693, -0.4143, 0.6981, \ + 0.9181, 2.9477, -0.3757, 0.8218, -0.5699, -0.5495, -0.0108, \ + 0.0854, 1.0189, 1.7904, 0.7577, -0.8434, -0.7300, -0.8675, \ + -0.6159, -0.3319, 1.4449, 0.6741}; + + struct s_ParamsAP Par; + + Par.Fs = NULL; Par.Fc = NULL; Par.Ns = NULL; Par.im = NULL; Par.ie = NULL; + + FILE *fid = NULL; + + + + /* Number of sample points */ + + n = 1024; + + + + /* Modulator (LP-random) */ + + m = (double*) calloc(n, sizeof(double)); + + if (m == NULL) {line=__LINE__-2; exitflag=-1; goto failed;} + + + for (i=0; i + +#include + +#include "f_ap_demodulation.c" + + + +#ifdef _WIN32 + + #define STR_NL "\r\n" + +#else + + #define STR_NL "\n" + +#endif + + + +int main(void) +{ + + /* Declarations and initializations */ + + int line; + + int exitflag = 0; + + + long i, j, iter, n; + + double dt, norm; + + double *m=NULL, *c=NULL, *s=NULL, *Ub=NULL, *t=NULL; + + double *out_m_wo_ub=NULL, *out_e_wo_ub=NULL, *out_E_wo_ub=NULL, \ + *out_m_w_ub=NULL, *out_e_w_ub=NULL, *out_E_w_ub=NULL; + + double w[] = {1.5648, 0.5312, 0.1413, 0.7588, -0.8616, -0.3586, 0.9106, -0.1787,\ + -0.0108, -0.0989, -0.3559, -0.4015, 0.2917, -0.3458, -1.1990, \ + 0.7651, -0.9884, -1.1668, 0.6584, -1.3693}; + + long ci[] = {0, 38, 46, 27, 35, 30, 29, 32, 36, 37, 41, 38, 45, 32, 50, 27, 45, \ + 38, 42, 30, 32, 48, 53, 35, 45, 50, 51}; + + struct s_ParamsAP Par; + + Par.Fs = NULL; Par.Fc = NULL; Par.Ns = NULL; Par.im = NULL; Par.ie = NULL; + + FILE *fid = NULL; + + + + /* Number of sample points */ + + n = 1024; + + + + /* Time step */ + + dt = 10/(double)(n-1); + + + + /* Modulator (LP-random) */ + + m = (double*) calloc(n, sizeof(double)); + + if (m == NULL) {line=__LINE__-2; exitflag=-1; goto failed;} + + + for (i=0; i + +#include + +#include "f_ap_demodulation.c" + + + +#ifdef _WIN32 + + #define STR_NL "\r\n" + +#else + + #define STR_NL "\n" + +#endif + + + +int main(void) +{ + + /* Declarations and initializations */ + + int line; + + int exitflag = 0; + + + long i, j, iter, n; + + double dt, smin, smmin; + + double *m1=NULL, *m2=NULL, *c1=NULL, *c2=NULL, *s=NULL, *Ub=NULL, *t=NULL; + + double *out_m1=NULL, *out_m2=NULL, *out_e1=NULL, *out_e2=NULL; + + double w[] = {1.5648, 0.5312, 0.1413, 0.7588, -0.8616, -0.3586, 0.9106, -0.1787,\ + -0.0108, -0.0989, -0.3559, -0.4015, 0.2917, -0.3458, -1.1990, \ + 0.7651, -0.9884, -1.1668, 0.6584, -1.3693, 0.7608, 0.7810, 0.9041,\ + 0.2338, 0.1767, 0.3911, 0.3206, 0.8155, 0.6135, 0.7600}; + + struct s_ParamsAP Par; + + Par.Fs = NULL; Par.Fc = NULL; Par.Ns = NULL; Par.im = NULL; Par.ie = NULL; + + FILE *fid = NULL; + + + + /* Number of sample points */ + + n = 1024; + + + + /* Time step */ + + dt = (double)10/(n-1); + + + + /* Modulators (LP-random) */ + + m1 = (double*) calloc(n, sizeof(double)); + + if (m1 == NULL) {line=__LINE__-2; exitflag=-1; goto failed;} + + + for (i=0; i=0; j--) + + m2[i] = m2[i] + w[2*j] * cos(2*M_PI*j*(n-i)/(double)n+w[1+2*j]); + + m2[i] = m2[i] + 2.581418146550079; + + m2[i] = m2[i] / 8.690964954126397 / 2; + } + + + + /* Carrier (regular spikes) */ + + c1 = (double*) calloc(n, sizeof(double)); + + if (c1 == NULL) {line=__LINE__-2; exitflag=-1; goto failed;} + + + c2 = (double*) calloc(n, sizeof(double)); + + if (c2 == NULL) {line=__LINE__-2; exitflag=-1; goto failed;} + + + for (i=4-1; i 1, demodulation is performed by + * using signal compression (see equations (12)-(13) in the paper cited + * in the header of this file. + * + * .Br - indicator of premature termination of the 'AP-Accelerated' algorithm + * when the λ factor drops below one. If .Br=1 (recommended), premature + * termination is assumed. Otherwise, if .Br=0, the AP-A is not stopped + * prematurely even when λ decreases below 1. This field is required + * only if .Al='A'. + * + * .im - array with the iteration numbers at which the modulator estimates + * have to be saved for the output. the first element is the length of + * the array (excluding the first element itself). + * + * .ie - array with the iteration numbers at which the infeasibility error, + * ϵ, values have to be saved for the output. The first element of .ie + * is the length of the array (excluding the first element itself). + * + * Two additional fields, .ns (number of sample points of the original + * signal) and .Nx (dimensions of the actual, possibly interpolated signal) + * are assigned values in this function. + * + * [Ub] - upper bound on the modulator. This array must have the same number of + * elements as the input signal or must be set to NULL (if no upper bound on + * the modulator is assumed). + * + * [t] - sampling coordinates of the input signal. This is a 2D array with the number + * of columns equal to the dimension of the input signal and the number of rows + * equal to the number of signal sample points. This input argument is optional + * and should be used only if the input signal is sampled + * nonuniformly. Otherwise, t must be set to NULL. If t is not NULL, then the + * input argument Par must contain the field .Nr (see above). + */ + +/* O U T P U T A R G U M E N T S + * + * [exitflag] - exit flag (non zero value indicates an error). + * + * [out_m#] - array with modulator estimates at algorithm iterations indicated by + * Par.im (memory allocated externally). The modulator estimates are + * arranged columnwise. + * + * [out_e#] - array with error estimates at algorithm iterations indicated by Par.ie + * (memory allocated externally). + * + * [iter#] - number of AP iterations. + */ + +/* N O N S T A N D A R D F U N C T I O N S U S E D + * + * (1) f_ap_compression, (2) f_ap_interpolation, (3) f_ap_s_Ub_init, + * + * (4) f_ap_mkl_dft_init, (5) l_ap_input_validation, (6) DftiFreeDescriptor, + * + * (7) f_ap_basic, (8) f_ap_accelerated, (9) f_ap_projected. + */ + + + +/***********************************************************************************/ +/*********************** DECLARATIONS AND INITIALIZATIONS *************************/ +/***********************************************************************************/ + + + int line; + + int exitflag = 0; + + + + /* Validation of the input data */ + + exitflag = f_ap_input_validation (s, Par, Ub, t); + + if (exitflag != 0) + + return exitflag; + + + + /* Basic variables */ + + long i; + + long nx = 1; + + long nx_2; + + long *ix_map = NULL; + + + DFTI_DESCRIPTOR_HANDLE dft_handle = 0; + + + double *s_local = NULL; + + double *Ub_local = NULL; + + double *s_local2 = NULL; + + double *Ub_local2 = NULL; + + double *s_local3 = NULL; + + double *Ub_local3 = NULL; + + double *pr_s = s; + + double *pr_Ub = Ub; + + + + /* Dimensions of the actual signal to be demodulated */ + + if (t != NULL) + + Par->Nx = Par->Nr; + + else + + Par->Nx = Par->Ns; + + + + /* Numbers of sample points */ + + if (t != NULL) + { + Par->ns = Par->Ns[0]; /* original signal */ + + for (i=0; i<(Par->D); i++) + + nx = nx * Par->Nx[i]; /* actual signal */ + } + + else + { + Par->ns = 1; + + for (i=0; i<(Par->D); i++) + { + Par->ns = Par->ns * Par->Ns[i]; /* original signal */ + + nx = nx * Par->Nx[i]; /* actual signal */ + } + } + + nx_2 = (nx / Par->Nx[Par->D-1]) * (Par->Nx[Par->D-1]+2-(Par->Nx[Par->D-1]%2)); + + + +/***********************************************************************************/ +/********************************* CALCULATION *************************************/ +/***********************************************************************************/ + + + /* Compression */ + + if (Par->Cp > 1) + { + + s_local = (double*) malloc(Par->ns*sizeof(double)); + + if (s_local == NULL) {line=__LINE__-2; exitflag=-1; goto failed;} + + + f_ap_compression (s_local, Par->ns, 1/(Par->Cp)); + + + if (Ub != NULL) + { + Ub_local = (double*) malloc(Par->ns*sizeof(double)); + + if (Ub_local == NULL) {line=__LINE__-2; exitflag=-1; goto failed;} + + + f_ap_compression (pr_Ub, Par->ns, 1/(Par->Cp)); + } + + else + + Ub_local = NULL; + + + pr_s = s_local; + + pr_Ub = Ub_local; + } + + + + /* Interpolation */ + + ix_map = (long*) malloc(Par->ns*sizeof(long)); + + if (ix_map == NULL) {line=__LINE__-2; exitflag=-1; goto failed;} + + + if (t != NULL) + { + s_local2 = (double*) malloc(nx*sizeof(double)); + + if (s_local2 == NULL) {line=__LINE__-2; exitflag=-1; goto failed;} + + + if (Ub != NULL) + { + Ub_local2 = (double*) malloc(nx*sizeof(double)); + + if (Ub_local2 == NULL) {line=__LINE__-2; exitflag=-1; goto failed;} + } + + else + + Ub_local2 = NULL; + + + exitflag = f_ap_interpolation (pr_s, Par, pr_Ub, t, s_local2, Ub_local2, \ + ix_map); + + if (exitflag != 0) goto finish; + + + pr_s = s_local2; + + pr_Ub = Ub_local2; + + } + + else + { + for (i=0; i<(Par->ns); i++) + + ix_map[i] = i; + } + + + + /* Local copies of the signal and upper bound arrays with the element placement + * replaced from the standard to the Intel MKL DFT's. The last dimension of the + * arrays is complemented by one or two elements as required by the MKL DFT + * routine. The mapping between the two conventions is provided as well. */ + + s_local3 = (double*) malloc(nx_2*sizeof(double)); + + if (s_local3 == NULL) {line=__LINE__-2; exitflag=-1; goto failed;} + + + if (Ub != NULL) + { + Ub_local3 = (double*) malloc(nx_2*sizeof(double)); + + if (Ub_local3 == NULL) {line=__LINE__-2; exitflag=-1; goto failed;} + } + else + + Ub_local3 = NULL; + + + exitflag = f_ap_s_Ub_init (pr_s, pr_Ub, ix_map, Par->D, Par->Nx, Par->ns, \ + s_local3, Ub_local3); + + if (exitflag != 0) goto finish; + + + pr_s = s_local3; + + pr_Ub = Ub_local3; + + + + /* Intel MKL DFT's descriptor */ + + exitflag = f_ap_mkl_dft_init (Par->D, Par->Nx, &dft_handle); + + if (exitflag != 0) goto finish; + + + + /* Demodulation */ + + if (Par->Al == 'B') + + exitflag = f_ap_basic (pr_s, Par, pr_Ub, ix_map, &dft_handle, out_m, out_e, \ + iter); + + else if (Par->Al == 'A') + + exitflag = f_ap_accelerated (pr_s, Par, pr_Ub, ix_map, &dft_handle, out_m, \ + out_e, iter); + + else if (Par->Al == 'P') + + exitflag = f_ap_projected (pr_s, Par, pr_Ub, ix_map, &dft_handle, out_m, \ + out_e, iter); + + if (exitflag != 0) goto finish; + + + + /* Decompression */ + + if (Par->Cp > 1) + { + f_ap_compression (out_m, Par->ns*(Par->im[0]), Par->Cp); + + pr_s = s_local; + } + + + +/***********************************************************************************/ +/************************** OUTPUT & MEMORY DEALLOCATION ***************************/ +/***********************************************************************************/ + + + goto finish; + + finish: + + free(s_local); + + free(Ub_local); + + free(ix_map); + + free(s_local2); + + free(Ub_local2); + + free(s_local3); + + free(Ub_local3); + + DftiFreeDescriptor (&dft_handle); + + return exitflag; + + failed: + + if (exitflag == -1) + + f_ap_print_error (__FUNCTION__, line, AP_ERR_MSG_MEM); + + goto finish; + + +} + + +#endif + diff --git a/C/library/h_ap_demodulation.h b/C/library/h_ap_demodulation.h new file mode 100644 index 0000000..90175ed --- /dev/null +++ b/C/library/h_ap_demodulation.h @@ -0,0 +1,196 @@ + +/* C O P Y R I G H T N O T I C E + * + * Copyright ©2021. Institute of Science and Technology Austria (IST Austria). + * All Rights Reserved. The underlying technology is protected by PCT Patent + * Application No. PCT/EP2021/054650. + * + * This file is part of the AP demodulation library, which is free software: you can + * redistribute it and/or modify it under the terms of the GNU General Public License + * as published by the Free Software Foundation in version 2. + * + * This program 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 v2 for more details. You + * should have received a copy of the GNU General Public License v2 along with this + * program. If not, see https://www.gnu.org/licenses/. + * + * Contact the Technology Transfer Office, IST Austria, Am Campus 1, + * A-3400 Klosterneuburg, Austria, +43-(0)2243 9000, twist@ist.ac.at, for commercial + * licensing opportunities. + * + * All other inquiries should be directed to the author, Mantas Gabrielaitis, + * mgabriel@ist.ac.at + */ + + +/* C O N T E N T S + * + * This is the header file for the f_ap_demodulation routine. It includes all needed + * external libraries, declares the input parameter structure for the + * f_ap_demodulation function, declares the prototype of the f_ap_demodulation + * function itself, defines macros with error mesages, and operating-system-dependent + * new-line control characters. + */ + + + +#ifndef ap_demodulation_h + + #define ap_demodulation_h + + + + /* EXTERNAL LIBRARIES */ + + #include + + #include + + #include + + #include + + + #include "mkl.h" + + #include "mkl_dfti.h" + + + + + /* INPUT PARAMETER STRUCTURE */ + + struct s_ParamsAP { + + double* Fs; + + double* Fc; + + char Al; + + double Et; + + long Ni; + + int D; + + long* Ns; + + long ns; + + long* Nr; + + long* Nx; + + double Cp; + + int Br; + + long* im; + + long* ie; + + }; + + + + + /* AP DEMODULATION FUNCTION */ + + int f_ap_demodulation (double*, struct s_ParamsAP*, double*, const double*, \ + double*, double*, long*); + + + + + /* ERROR MESSAGES (INVALID INPUT)*/ + + #define AP_ERR_MSG_AL "AP algorithm, set by Par.Al, must be either \'B\' "\ + "(Basic), \'A\' (Accelerated), or \'P\' (Projected)!" + + #define AP_ERR_MSG_D "The number of signal dimensions, set by Par.D, must be "\ + "equal to 1, 2, or 3!" + + #define AP_ERR_MSG_FS "Sampling frequencies, set by Par.Fs, must be positive "\ + "real numbers!" + + #define AP_ERR_MSG_FC "Cutoff frequencies, set by Par.Fc, must be "\ + "non-negative real numbers!" + + #define AP_ERR_MSG_FC2 "Cutoff frequencies, set by Par.Fc, cannot be higher "\ + "than half of the sampling frequencies set by Par.Fs!" + + #define AP_ERR_MSG_ET "Error tolerance, set by Par.Et, must be a real number!" + + #define AP_ERR_MSG_NI "Iteration number, set by Par.Niter, must be a positive "\ + "number!" + + #define AP_ERR_MSG_NS "Numbers of elements of the signal array in each "\ + "dimension, set by Par.Ns, must be integers higher than 1!" + + #define AP_ERR_MSG_NR "Numbers of elements of the interpolated signal array in "\ + "each dimension, set by Par.Nr, must be integers higher "\ + "than 1!" + + #define AP_ERR_MSG_CP "Compression parameter, set by Par.Cp, must be a real "\ + "number not smaller than 1!" + + #define AP_ERR_MSG_BR "Indicator of premature termination of the AP-A "\ + "algorithm, set by Par.Br, must be either 0 or 1!" + + #define AP_ERR_MSG_IE "The number of error-sampling iterations, set by "\ + "Par.ie[0], must be positive!" + + #define AP_ERR_MSG_IE2 "Error-sampling iterations, set by Par.ie, must be "\ + "nonnegative strictly increasing numbers!" + + #define AP_ERR_MSG_IM "The number of modulator-sampling iterations, set by "\ + "Par.im[0], must be positive!" + + #define AP_ERR_MSG_IM2 "Modulator-sampling iterations, set by Par.im, must be "\ + "nonnegative strictly increasing numbers!" + + #define AP_ERR_MSG_S "Array with the signal values, set by the input argument "\ + "s, must consist of real numbers!" + + #define AP_ERR_MSG_UB "Array with the upper bound values for modulator, set by "\ + "the input argument Ub, must consist of real numbers not "\ + "smaller than entries of the absolute value signal!" + + #define AP_ERR_MSG_T "Array with the sampling coordinates, set by the input "\ + "argument t, must consist of real numbers!" + + + + + /* ERROR MESSAGES (OTHER)*/ + + #define AP_ERR_MSG_PRE "(!) ERROR in AP Demodulation" + + #define AP_ERR_MSG_MEM "Out of memmory!" + + + + + + /* NEW-LINE CONTROL CHARACTERS */ + + #ifdef _WIN32 + + #define AP_STR_NL "\r\n" + + #define AP_STR_NL2 "\r\n\r\n" + + #else + + #define AP_STR_NL "\n" + + #define AP_STR_NL2 "\n\n" + + #endif + + +#endif + + \ No newline at end of file diff --git a/C/library/l_ap_algorithms.c b/C/library/l_ap_algorithms.c new file mode 100644 index 0000000..c1b9c43 --- /dev/null +++ b/C/library/l_ap_algorithms.c @@ -0,0 +1,1125 @@ + +/* C O P Y R I G H T N O T I C E + * + * Copyright ©2021. Institute of Science and Technology Austria (IST Austria). + * All Rights Reserved. The underlying technology is protected by PCT Patent + * Application No. PCT/EP2021/054650. + * + * This file is part of the AP demodulation library, which is free software: you can + * redistribute it and/or modify it under the terms of the GNU General Public License + * as published by the Free Software Foundation in version 2. + * + * This program 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 v2 for more details. You + * should have received a copy of the GNU General Public License v2 along with this + * program. If not, see https://www.gnu.org/licenses/. + * + * Contact the Technology Transfer Office, IST Austria, Am Campus 1, + * A-3400 Klosterneuburg, Austria, +43-(0)2243 9000, twist@ist.ac.at, for commercial + * licensing opportunities. + * + * All other inquiries should be directed to the author, Mantas Gabrielaitis, + * mgabriel@ist.ac.at + */ + + +/* C O N T E N T S + * + * Three functions implementing different alternating projection algorithms of + * amplitude demodulation: + * + * (1) f_ap_basic, + * + * (2) f_ap_accelerated, + * + * (3) f_ap_projected. + */ + + + +int f_ap_basic (double* s, const struct s_ParamsAP* Par, const double* Ub, \ + const long* ix_map, DFTI_DESCRIPTOR_HANDLE* dft_handle, double* m_out, \ + double* e_out, long* iter) +{ +/* S H O R T D E S C R I P T I O N + * + * Calculates the modulator of a signal by using the AP-Basic algorithm. Signals + * defined in up to 3 dimensions are allowed. + */ + +/* I N P U T A R G U M E N T S + * + * [s] - input signal + 2 additional array elements along the last dimension. This + * input argument is modified in-place! + * + * [Par] - pointer to the structure with demodulation parameters: + * + * .D - number of signal dimensions. + * + * .Fs - sampling frequencies for each dimension of the signal. This is an + * array of D elements. + * + * .Fc - cutoff frequencies of the modulator for each dimension of the + * signal. This is an array of D elements. + * + * .Et - infeasibility error tolerance used to control the termination of the + * demodulation algorithm. The iterative process is stopped when the + * infeasibility error, ϵ, drops to the level of .Et or below. If + * .Et ≤ 0, then the maximum allowed number of iterations, .Ni (see + * below), is completed. + * + * .Ni - maximum number of allowed iterations of the chosen AP algorithm. The + * chosen AP algorithm is iterated not more than .Ni times independent + * of whether ϵ drops at or below .Et. + * + * .ns - number of sample points of the original input signal. + * + * .Nx - number of elements of the provided input signal (possibly + * interpolated) in every dimension. This is an array with the number + * of elements equal to the number of signal dimensions. It does not + * assume the additional two elements in the last dimension of s. + * + * .im - array with the iteration numbers at which the modulator estimates + * have to be saved for the output. the first element is the length of + * the array (excluding the first element itself). + * + * .ie - array with the iteration numbers at which the infeasibility error, + * ϵ, values have to be saved for the output. The first element of .ie + * is the length of the array (excluding the first element itself). + * + * [Ub] - upper bound on the modulator. This array must have the same number of + * elements as the input signal (does not include the additional two elements + * in the last dimension of s). + * + * [ix_map] - indexes of the moddulator elements to be saved for the output. This + * array is either NULL or consists of the same number of elements as the + * original input signal (before any possible interpolation). + * + * [dft_handle] - address of the comitted descriptor handle of the Intel MKL DFT. + */ + +/* O U T P U T A R G U M E N T S + * + * [exitflag] - exit flag. + * + * [m_out#] - array with modulator estimates (memory allocated externally). + * + * [e_out#] - array with error estimates (memory allocated externally). + * + * [iter#] - the actual number of iterations used. + */ + +/* N O N S T A N D A R D F U N C T I O N S U S E D + * + * (1) f_ap_abs_scaled_max_abs, (2) f_ap_mkl_dft_PMw. + */ + + +/***********************************************************************************/ +/*********************** DECLARATIONS AND INITIALIZATIONS *************************/ +/***********************************************************************************/ + + int line; + + int exitflag = 0; + + + long i; + + long i_aux; + + long iter_m = 1; + + long iter_e = 1; + + long nx = 1; + + long nx_2; + + long *iL = NULL; + + long *iR = NULL; + + + double E; + + double Etol; + + double max_s_abs; + + double s_old; + + double aux; + + + double *s_abs = NULL; + + + + /* Number of sample points of the provided signal */ + + for (i=0; i<(Par->D); i++) + + nx = nx * Par->Nx[i]; + + nx_2 = (nx / Par->Nx[Par->D-1]) * (Par->Nx[Par->D-1]+2-(Par->Nx[Par->D-1]%2)); + + + + /* Indexes of the left and right cutoff frequencies */ + + iL = (long*) malloc((Par->D)*sizeof(long)); + + if (iL == NULL) {line=__LINE__-2; exitflag=-1; goto failed;} + + + iR = (long*) malloc((Par->D)*sizeof(long)); + + if (iR == NULL) {line=__LINE__-2; exitflag=-1; goto failed;} + + + for (i=0; i<(Par->D); i++) + { + iL[i] = 1 + (long) ceil(Par->Fc[i] / (Par->Fs[i] / Par->Nx[i])); + + iR[i] = Par->Nx[i] - iL[i]; + } + + + + /* Normalized absolute-value version of the signal */ + + s_abs = (double*) malloc(nx_2*sizeof(double)); + + if (s_abs == NULL) {line=__LINE__-2; exitflag=-1; goto failed;} + + + max_s_abs = f_ap_abs_scaled_max_abs (s, nx_2, s_abs); + + + + /* Initialization of the error tolerance variable */ + + if (Par->Et > 0) + + Etol = (Par->Et / max_s_abs) * (Par->Et / max_s_abs) * nx; + + else + + Etol = Par->Et; + + + + /* Initialization of the modulator and infeasibility error variables */ + + E = 0; + + for (i=0; iim[iter_m] == 0) + { + for (i=0; i<(Par->ns); i++) + + m_out[i] = s_abs[ix_map[i]] * max_s_abs; + + iter_m = iter_m + 1; + } + + + + /* Readout of the infeasibility error */ + + if (Par->ie[iter_e] == 0) + { + e_out[0] = sqrt(E / nx); + + iter_e = iter_e + 1; + } + + + +/***********************************************************************************/ +/********************************* CALCULATION *************************************/ +/***********************************************************************************/ + + + /* Alternating projections */ + + *iter = 0; + + while (E > Etol && Par->Ni > *iter) + { + *iter = *iter + 1; + + + + /* Projection onto the set Mw */ + + exitflag = f_ap_mkl_dft_PMw (s, Par->D, Par->Nx, iL, iR, dft_handle); + + if (exitflag != 0) goto finish; + + + + /* Projection onto the set Cd; error estimate */ + + E = 0; + + for (i=0; i Ub[i]) + + s[i] = Ub[i]; + + + + /* Infeasibility error */ + + aux = (s[i]-s_old); + + E = E + aux * aux; + + } + + + + /* Output (modulator) */ + + if ( iter_m <= Par->im[0] && (*iter == Par->im[iter_m] || \ + (E <= Etol && Par->im[0] == 1 && Par->im[1] == Par->Ni)) ) + { + i_aux = (iter_m-1)*(Par->ns); + + for (i=0; i<(Par->ns); i++) + + m_out[i+i_aux] = s[ix_map[i]] * max_s_abs; + + iter_m = iter_m + 1; + } + + + + /* Output (infeasibility error) */ + + if ( iter_e <= Par->ie[0] && (*iter == Par->ie[iter_e] || \ + (E <= Etol && Par->ie[0] == 1 && Par->ie[1] == Par->Ni)) ) + { + e_out[iter_e-1] = sqrt(E / nx); + + iter_e = iter_e + 1; + } + + } + + + +/***********************************************************************************/ +/************************** OUTPUT & MEMORY DEALLOCATION ***************************/ +/***********************************************************************************/ + + + goto finish; + + finish: + + free(iL); + + free(iR); + + free(s_abs); + + return exitflag; + + failed: + + if (exitflag == -1) + + f_ap_print_error (__FUNCTION__, line, AP_ERR_MSG_MEM); + + goto finish; + + +} + + + + +int f_ap_accelerated (double* s, const struct s_ParamsAP* Par, const double* Ub, \ + const long* ix_map, DFTI_DESCRIPTOR_HANDLE* dft_handle, double* m_out, \ + double* e_out, long* iter) +{ +/* S H O R T D E S C R I P T I O N + * + * Calculates the modulator of a signal by using the AP-Accelerated algorithm. + * Signals defined in up to 3 dimensions are allowed. + */ + +/* I N P U T A R G U M E N T S + * + * [s] - input signal + 2 additional array elements along the last dimension. This + * input argument is modified in-place! + * + * [Par] - pointer to the structure with demodulation parameters: + * + * .Fs - sampling frequencies. This is an array of D elements with sampling + * frequency for each dimension of the signal. + * + * .Fc - cutoff frequencies of the modulator. This is an array of D elements + * with cutoff frequencies for each dimension of the signal. + * + * .Et - infeasibility error tolerance used to control the termination of the + * demodulation algorithm. The iterative process is stopped when the + * infeasibility error, ϵ, drops to the level of .Et or below. If + * .Et ≤ 0, then the maximum allowed number of iterations, .Ni (see + * below), is completed. + * + * .Ni - maximum number of allowed iterations of the chosen AP algorithm. The + * chosen AP algorithm is iterated not more than .Ni times independent + * of whether ϵ drops at or below .Et. + * + * .D - number of signal dimensions. + * + * .ns - number of sample points of the original input signal. + * + * .Nx - number of elements of the provided input signal (possibly + * interpolated) in every dimension. This is an array with the number + * of elements equal to the number of signal dimensions. It does not + * assume the additional two elements in the last dimension of s. + * + * .Br - indicator of premature termination of the 'AP-Accelerated' algorithm + * when the λ factor drops below one. If .Br=1 (recommended), premature + * termination is assumed. Otherwise, if .Br=0, the AP-A is not stopped + * prematurely even when λ decreases below 1. + * + * .im - array with the iteration numbers at which the modulator estimates + * have to be saved for the output. the first element is the length of + * the array (excluding the first element itself). + * + * .ie - array with the iteration numbers at which the infeasibility error, + * ϵ, values have to be saved for the output. The first element of .ie + * is the length of the array (excluding the first element itself). + * + * [Ub] - upper bound on the modulator. This array must have the same number of + * elements as the input signal (does not include the additional two elements + * in the last dimension of s). + * + * [ix_map] - indexes of the moddulator elements to be saved for the output. This + * array is either NULL or consists of the same number of elements as the + * original input signal (before any possible interpolation). + * + * [dft_handle] - address of the comitted descriptor handle of the Intel MKL DFT. + */ + +/* O U T P U T A R G U M E N T S + * + * [exitflag] - exit flag. + * + * [m_out#] - array with modulator estimates (memory allocated externally). + * + * [e_out#] - array with error estimates (memory allocated externally). + * + * [iter#] - the actual number of iterations used. + */ + +/* N O N S T A N D A R D F U N C T I O N S U S E D + * + * (1) f_ap_abs_scaled_max_abs, (2) f_ap_mkl_dft_PMw. + */ + + +/***********************************************************************************/ +/*********************** DECLARATIONS AND INITIALIZATIONS *************************/ +/***********************************************************************************/ + + + int line; + + int exitflag = 0; + + + long i; + + long i_aux; + + long iter_m = 1; + + long iter_e = 1; + + long nx = 1; + + long nx_2; + + long *iL = NULL; + + long *iR = NULL; + + + double E; + + double Etol; + + double max_s_abs; + + double lambda; + + double nom; + + double denom; + + + double *a = NULL; + + double *b = NULL; + + double *s_abs = NULL; + + + + /* Number of sample points of the provided signal */ + + for (i=0; i<(Par->D); i++) + + nx = nx * Par->Nx[i]; + + nx_2 = (nx / Par->Nx[Par->D-1]) * (Par->Nx[Par->D-1]+2-(Par->Nx[Par->D-1]%2)); + + + + /* Indexes of the left and right cutoff frequencies */ + + iL = (long*) malloc((Par->D)*sizeof(long)); + + if (iL == NULL) {line=__LINE__-2; exitflag=-1; goto failed;} + + + iR = (long*) malloc((Par->D)*sizeof(long)); + + if (iR == NULL) {line=__LINE__-2; exitflag=-1; goto failed;} + + + for (i=0; i<(Par->D); i++) + { + iL[i] = 1 + (long) ceil(Par->Fc[i] / (Par->Fs[i] / Par->Nx[i])); + + iR[i] = Par->Nx[i] - iL[i]; + } + + + + /* Normalized absolute-value version of the signal */ + + s_abs = (double*) malloc(nx_2*sizeof(double)); + + if (s_abs == NULL) {line=__LINE__-2; exitflag=-1; goto failed;} + + + max_s_abs = f_ap_abs_scaled_max_abs (s, nx_2, s_abs); + + + + /* Initialization of the error tolerance variable */ + + if (Par->Et > 0) + + Etol = (Par->Et / max_s_abs) * (Par->Et / max_s_abs) * nx; + + else + + Etol = Par->Et; + + + + /* Initialization of the modulator-related and infeasibility error variables */ + + a = (double*) calloc(nx_2,sizeof(double)); + + if (a == NULL) {line=__LINE__-2; exitflag=-1; goto failed;} + + + b = (double*) malloc(nx_2*sizeof(double)); + + if (b == NULL) {line=__LINE__-2; exitflag=-1; goto failed;} + + + nom = 0; + + E = 0; + + for (i=0; iim[iter_m] == 0) + { + for (i=0; i<(Par->ns); i++) + + m_out[i] = s_abs[ix_map[i]] * max_s_abs; + + iter_m = iter_m + 1; + } + + + + /* Readout of the infeasibility error */ + + if (Par->ie[iter_e] == 0) + { + e_out[0] = sqrt(E / nx); + + iter_e = iter_e + 1; + } + + + +/***********************************************************************************/ +/********************************* CALCULATION *************************************/ +/***********************************************************************************/ + + + /* Alternating projections */ + + *iter = 0; + + while (E > Etol && Par->Ni > *iter) + { + *iter = *iter + 1; + + + + /* Projection onto the set Mw */ + + exitflag = f_ap_mkl_dft_PMw (b, Par->D, Par->Nx, iL, iR, dft_handle); + + if (exitflag != 0) goto finish; + + + + /* Factor lambda */ + + denom = 0; + + for (i=0; iBr != 0) + + break; + + + + /* Projection onto the set Cd; a, b, nom, and error estimates */ + + nom = 0; + + E = 0; + + for (i=0; i Ub[i]) + + s[i] = Ub[i]; + + + + /* Nominator of the factor lambda for the next iteration */ + + b[i] = s[i] - a[i]; + + nom = nom + b[i] * b[i]; + + + + /* Infeasibility error */ + + E = nom; + } + + + + /* Output (modulator) */ + + if ( iter_m <= Par->im[0] && (*iter == Par->im[iter_m] || \ + (E <= Etol && Par->im[0] == 1 && Par->im[1] == Par->Ni)) ) + { + i_aux = (iter_m-1)*(Par->ns); + + for (i=0; i<(Par->ns); i++) + + m_out[i+i_aux] = s[ix_map[i]] * max_s_abs; + + iter_m = iter_m + 1; + } + + + + /* Output (infeasibility error) */ + + if ( iter_e <= Par->ie[0] && (*iter == Par->ie[iter_e] || \ + (E <= Etol && Par->ie[0] == 1 && Par->ie[1] == Par->Ni)) ) + { + e_out[iter_e-1] = sqrt(E / nx); + + iter_e = iter_e + 1; + } + } + + + +/***********************************************************************************/ +/************************** OUTPUT & MEMORY DEALLOCATION ***************************/ +/***********************************************************************************/ + + + goto finish; + + finish: + + free(iL); + + free(iR); + + free(s_abs); + + free(a); + + free(b); + + return exitflag; + + failed: + + if (exitflag == -1) + + f_ap_print_error (__FUNCTION__, line, AP_ERR_MSG_MEM); + + goto finish; + + +} + + + + +int f_ap_projected (double* s, const struct s_ParamsAP* Par, const double* Ub, \ + const long* ix_map, DFTI_DESCRIPTOR_HANDLE* dft_handle, double* m_out, \ + double* e_out, long* iter) +{ +/* S H O R T D E S C R I P T I O N + * + * Calculates the modulator of a signal by using the AP-Basic algorithm. Signals + * defined in up to 3 dimensions are allowed. + */ + +/* I N P U T A R G U M E N T S + * + * [s] - input signal + 2 additional array elements along the last dimension. This + * input argument is modified in-place! + * + * [Par] - pointer to the structure with demodulation parameters: + * + * .Fs - sampling frequencies. This is an array of D elements with sampling + * frequency for each dimension of the signal. + * + * .Fc - cutoff frequencies of the modulator. This is an array of D elements + * with cutoff frequencies for each dimension of the signal. + * + * .Et - infeasibility error tolerance used to control the termination of the + * demodulation algorithm. The iterative process is stopped when the + * infeasibility error, ϵ, drops to the level of .Et or below. If + * .Et ≤ 0, then the maximum allowed number of iterations, .Ni (see + * below), is completed. + * + * .Ni - maximum number of allowed iterations of the chosen AP algorithm. The + * chosen AP algorithm is iterated not more than .Ni times independent + * of whether ϵ drops at or below .Et. + * + * .D - number of signal dimensions. + * + * .ns - number of sample points of the original input signal. + * + * .Nx - number of elements of the provided input signal (possibly + * interpolated) in every dimension. This is an array with the number + * of elements equal to the number of signal dimensions. It does not + * assume the additional two elements in the last dimension of s. + * + * .im - array with the iteration numbers at which the modulator estimates + * have to be saved for the output. the first element is the length of + * the array (excluding the first element itself). + * + * .ie - array with the iteration numbers at which the infeasibility error, + * ϵ, values have to be saved for the output. The first element of .ie + * is the length of the array (excluding the first element itself). + * + * [Ub] - upper bound on the modulator. This array must have the same number of + * elements as the input signal (does not include the additional two elements + * in the last dimension of s). + * + * [ix_map] - indexes of the moddulator elements to be saved for the output. This + * array is either NULL or consists of the same number of elements as the + * original input signal (before any possible interpolation). + * + * [dft_handle] - address of the comitted descriptor handle of the Intel MKL DFT. + */ + +/* O U T P U T A R G U M E N T S + * + * [exitflag] - exit flag. + * + * [m_out#] - array with modulator estimates (memory allocated externally). + * + * [e_out#] - array with error estimates (memory allocated externally). + * + * [iter#] - the actual number of iterations used. + */ + +/* N O N S T A N D A R D F U N C T I O N S U S E D + * + * (1) f_ap_abs_scaled_max_abs, (2) f_ap_mkl_dft_PMw. + */ + + +/***********************************************************************************/ +/*********************** DECLARATIONS AND INITIALIZATIONS *************************/ +/***********************************************************************************/ + + + int line; + + int exitflag = 0; + + + long i; + + long i_aux; + + long iter_m = 1; + + long iter_e = 1; + + long nx = 1; + + long nx_2; + + long *iL = NULL; + + long *iR = NULL; + + + double E; + + double Etol; + + double max_s_abs; + + double aux; + + double aux2; + + + double *s_abs = NULL; + + double *a = NULL; + + double *c = NULL; + + + + /* Number of sample points of the provided signal */ + + for (i=0; i<(Par->D); i++) + + nx = nx * Par->Nx[i]; + + nx_2 = (nx / Par->Nx[Par->D-1]) * (Par->Nx[Par->D-1]+2-(Par->Nx[Par->D-1]%2)); + + + + /* Indexes of the left and right cutoff frequencies */ + + iL = (long*) malloc((Par->D)*sizeof(long)); + + if (iL == NULL) {line=__LINE__-2; exitflag=-1; goto failed;} + + + iR = (long*) malloc((Par->D)*sizeof(long)); + + if (iR == NULL) {line=__LINE__-2; exitflag=-1; goto failed;} + + + for (i=0; i<(Par->D); i++) + { + iL[i] = 1 + (long) ceil(Par->Fc[i] / (Par->Fs[i] / Par->Nx[i])); + + iR[i] = Par->Nx[i] - iL[i]; + } + + + + /* Normalized absolute-value version of the signal */ + + s_abs = (double*) malloc(nx_2*sizeof(double)); + + if (s_abs == NULL) {line=__LINE__-2; exitflag=-1; goto failed;} + + + max_s_abs = f_ap_abs_scaled_max_abs (s, nx_2, s_abs); + + + + /* Initialization of the error tolerance variable */ + + if (Par->Et > 0) + + Etol = (Par->Et / max_s_abs) * (Par->Et / max_s_abs) * nx; + + else + + Etol = Par->Et; + + + + /* Initialization of the infeasibility error and modulator-related variables */ + + a = (double*) malloc(nx_2*sizeof(double)); + + if (a == NULL) {line=__LINE__-2; exitflag=-1; goto failed;} + + + c = (double*) malloc(nx_2*sizeof(double)); + + if (c == NULL) {line=__LINE__-2; exitflag=-1; goto failed;} + + + E = 0; + + for (i=0; iim[iter_m] == 0) + { + for (i=0; i<(Par->ns); i++) + + m_out[i] = s_abs[ix_map[i]] * max_s_abs; + + iter_m = iter_m + 1; + } + + + + /* Readout of the infeasibility error */ + + if (Par->ie[iter_e] == 0) + { + e_out[0] = sqrt(E / nx); + + iter_e = iter_e + 1; + } + + + +/***********************************************************************************/ +/********************************* CALCULATION *************************************/ +/***********************************************************************************/ + + + /* Alternating projections */ + + *iter = 0; + + while (E > Etol && Par->Ni > *iter) + { + *iter = *iter + 1; + + + /* Projection onto the set Mw */ + + exitflag = f_ap_mkl_dft_PMw (a, Par->D, Par->Nx, iL, iR, dft_handle); + + if (exitflag != 0) goto finish; + + + + /* Projection onto the set Cd; error estimate */ + + E = 0; + + for (i=0; i Ub[i]) + + s[i] = Ub[i]; + + + + /* c and a */ + + aux2 = s[i] - a[i]; + + c[i] = c[i] + aux2; + + a[i] = s[i]; + + + + /* Infeasibility error */ + + E = E + aux * aux + aux2 * aux2; + + } + + + + /* Output (modulator) */ + + if ( iter_m <= Par->im[0] && (*iter == Par->im[iter_m] || \ + (E <= Etol && Par->im[0] == 1 && Par->im[1] == Par->Ni)) ) + { + i_aux = (iter_m-1)*(Par->ns); + + for (i=0; i<(Par->ns); i++) + + m_out[i+i_aux] = s[ix_map[i]] * max_s_abs; + + iter_m = iter_m + 1; + } + + + + /* Output (infeasibility error) */ + + if ( iter_e <= Par->ie[0] && (*iter == Par->ie[iter_e] || \ + (E <= Etol && Par->ie[0] == 1 && Par->ie[1] == Par->Ni)) ) + { + e_out[iter_e-1] = sqrt(E / nx / 2); + + iter_e = iter_e + 1; + } + + } + + + +/***********************************************************************************/ +/************************** OUTPUT & MEMORY DEALLOCATION ***************************/ +/***********************************************************************************/ + + + goto finish; + + finish: + + free(iL); + + free(iR); + + free(s_abs); + + free(a); + + free(c); + + return exitflag; + + failed: + + if (exitflag == -1) + + f_ap_print_error (__FUNCTION__, line, AP_ERR_MSG_MEM); + + goto finish; + + +} + diff --git a/C/library/l_ap_auxiliary.c b/C/library/l_ap_auxiliary.c new file mode 100644 index 0000000..eedd100 --- /dev/null +++ b/C/library/l_ap_auxiliary.c @@ -0,0 +1,1231 @@ + +/* C O P Y R I G H T N O T I C E + * + * Copyright ©2021. Institute of Science and Technology Austria (IST Austria). + * All Rights Reserved. The underlying technology is protected by PCT Patent + * Application No. PCT/EP2021/054650. + * + * This file is part of the AP demodulation library, which is free software: you can + * redistribute it and/or modify it under the terms of the GNU General Public License + * as published by the Free Software Foundation in version 2. + * + * This program 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 v2 for more details. You + * should have received a copy of the GNU General Public License v2 along with this + * program. If not, see https://www.gnu.org/licenses/. + * + * Contact the Technology Transfer Office, IST Austria, Am Campus 1, + * A-3400 Klosterneuburg, Austria, +43-(0)2243 9000, twist@ist.ac.at, for commercial + * licensing opportunities. + * + * All other inquiries should be directed to the author, Mantas Gabrielaitis, + * mgabriel@ist.ac.at + */ + + +/* C O N T E N T S + * + * Eight auxiliary functions for amplitude demodulation via alternating projections: + * + * (1) f_ap_print_error, + * + * (2) f_ap_minmax, + * + * (3) f_ap_abs_scaled_max_abs, + * + * (4) f_ap_compression, + * + * (5) f_ap_interpolation, + * + * (6) f_ap_s_Ub_init, + * + * (7) f_ap_mkl_dft_init, + * + * (8) f_ap_mkl_dft_PMw. + */ + + + +void f_ap_print_error (const char* function, const int line, const char* message) +{ +/* S H O R T D E S C R I P T I O N + * + * Prints error message to stderr or MATLAB command window. + */ + +/* I N P U T A R G U M E N T S + * + * [function] - name of the function where the error is detected. + * + * [line] - number of line in the code of the dunction where error is detected. + * + * [message] - error message. + */ + +/* O U T P U T A R G U M E N T S + * + * None. + */ + + + #ifdef AP_MEX_COMPILATION + + mexPrintf (AP_STR_NL AP_ERR_MSG_PRE " (line %d in %s): " AP_STR_NL2 \ + "%s" AP_STR_NL2, line, function, message); + + #else + + fprintf (stderr, AP_STR_NL AP_ERR_MSG_PRE " (line %d in %s): " AP_STR_NL2 \ + "%s" AP_STR_NL2, line, function, message); + + #endif + +} + + + + +void f_ap_minmax (const double* in, const long N, double* out_min, double* out_max) +{ +/* S H O R T D E S C R I P T I O N + * + * Finds the minimum and maximum values of the input array. + */ + +/* I N P U T A R G U M E N T S + * + * [in] - input array. + * + * [N] - number of elements in the array of the input signal. + */ + +/* O U T P U T A R G U M E N T S + * + * [out_min#] - the minimum value of the input array (this is the address of an + * externally defined scalar variable. + * + * [out_max#] - the maximum value of the input array (this is the address of an + * externally defined scalar variable. + */ + + + /* Initializations and declarations */ + + long i; + + + out_min[0] = in[0]; + + out_max[0] = in[0]; + + + + /* Calculation */ + + for (i=1; i in[i]) + + out_min[0] = in[i]; + + else if (out_max[0] < in[i]) + + out_max[0] = in[i]; + } + +} + + + + +double f_ap_abs_scaled_max_abs (const double* in, const long N, double* out) +{ +/* S H O R T D E S C R I P T I O N + * + * Calculates the absolute value of the input signal normzlized by its maximum. + */ + +/* I N P U T A R G U M E N T S + * + * [in] - input signal. + * + * [N] - number of elements in the array of the input signal. + */ + +/* O U T P U T A R G U M E N T S + * + * [out#] - normalized absolute-value signal (memory allocated externally); + * + * [maxval] - absolute maximum value of the input signal. + */ + + + /* Initializations and declarations */ + + long i; + + double maxval = 0; + + + + /* The maximum absolute value of the signal */ + + for (i=0; i maxval) + + maxval = out[i]; + } + + + + /* Scaled signal */ + + for (i=0; i0)-(s[i]<0)) * pow(fabs(s[i]),p); + +} + + + + +int f_ap_interpolation (const double* s, const struct s_ParamsAP* Par, \ + const double* Ub, const double* t, double* s_out, double* Ub_out, \ + long* ix_out) +{ +/* S H O R T D E S C R I P T I O N + * + * Interpolates the input signal on a refined uniform grid following Eq. 23 in + * M. Gabrielaitis IEEE Trans. Signal Process., vol. 69, pp. 4039-4054, 2021. + */ + +/* I N P U T A R G U M E N T S + * + * [s] - input signal. + * + * [Par] - pointer to the structure with interpolation parameters: + * + * .D - number of the input signal dimensions. + * + * .ns - number of sample points of the original input signal. + * + * .Nr - number of sample points on the refined uniform grid. This is an + * array of D elements. + * + * [Ub] - array with the upper bound values on the modulator. This is an array of the + * same size as the input signal s or Ub == NULL. In the latter case, no upper + * bound is assumed. + * + * [t] - sampling coordinates of the input signal. This is a 2D array with the number + * of columns equal to the dimension of the input signal (Par.D) and the number + * of rows equal to the total number of signal sample points, Par.Ns[0]. + */ + +/* O U T P U T A R G U M E N T S + * + * [exitflag] - exit flag. + * + * [s_out#] - interpolated signal. The number of dimensions of this array is given by + * Par.D. + * + * [Ub_out#] - interpolated upper bounds on the modulator. The dimensions of this + * array are the same as of s_out. + * + * [ix_out#] - linear indexes of the elements of s_out corresponding to every' + * element of s. + */ + +/* N O N S T A N D A R D F U N C T I O N S U S E D + * + * (1) f_ap_print_error, (2) f_ap_minmax. + */ + + + /* Initializations and declarations */ + + int line; + + int exitflag = 0; + + + int log_aux; + + long i1, i2, i3; + + long nr = 1; + + long *cumnr = NULL; + + long ix; + + long *ix_aux = NULL; + + + double *tmin = NULL; + + double *tmax = NULL; + + double *dt = NULL; + + double r2; + + double *r2_all = NULL; + + + + + /* Memory allocation for some of the internal arrays */ + + cumnr = (long*) malloc((Par->D)*sizeof(long)); + + if (cumnr == NULL) {line=__LINE__-2; exitflag=-1; goto failed;} + + + ix_aux = (long*) malloc((Par->D)*sizeof(long)); + + if (ix_aux == NULL) {line=__LINE__-2; exitflag=-1; goto failed;} + + + tmin = (double*) malloc((Par->D)*sizeof(double)); + + if (tmin == NULL) {line=__LINE__-2; exitflag=-1; goto failed;} + + + tmax = (double*) malloc((Par->D)*sizeof(double)); + + if (tmax == NULL) {line=__LINE__-2; exitflag=-1; goto failed;} + + + dt = (double*) malloc((Par->D)*sizeof(double)); + + if (dt == NULL) {line=__LINE__-2; exitflag=-1; goto failed;} + + + + /* Number of sample points in the original and interpolated arrays */ + + for (i1=0; i1<(Par->D); i1++) + { + nr = nr * Par->Nr[i1]; + + if (i1==0) + + cumnr[i1] = Par->Nr[i1]; + + else + + cumnr[i1] = Par->Nr[i1] * cumnr[i1-1]; + } + + + + /* Memory allocation for r2_all arrays */ + + r2_all = (double*) malloc((Par->ns)*sizeof(double)); + + if (r2_all == NULL) {line=__LINE__-2; exitflag=-1; goto failed;} + + + + /* Preparation of the signal and upper bound output arrays */ + + for (i1=0; i1D); i1++) + { + f_ap_minmax (t+(i1*(Par->ns)), (Par->ns), tmin+i1, tmax+i1); + + dt[i1] = (tmax[i1] - tmin[i1]) / (Par->Nr[i1] - 1); + } + + + + /* Interpolation */ + + for (i1=0; i1<(Par->ns); i1++) + { + /* The closest point on the new grid and the corresponding distance */ + + r2 = 0; + + for (i2=0; i2<(Par->D); i2++) + { + i3 = i1 + i2*(Par->ns); + + ix_aux[i2] = lround((t[i3]-tmin[i2]) / dt[i2]); + + r2 = r2 + (t[i3]-tmin[i2]-ix_aux[i2]*dt[i2]) * \ + (t[i3]-tmin[i2]-ix_aux[i2]*dt[i2]); + } + + ix = ix_aux[0]; + + for (i2=1; i2<(Par->D); i2++) + + ix = ix + ix_aux[i2] * cumnr[i2-1]; + + + /* Checking whether the closest point on the new grid has already been used*/ + + log_aux = 0; + + for (i2=0; i2ns); i1++) + + ix_out[i1] = labs(ix_out[i1]); + + + + /* Output & Memory deallocation */ + + goto finish; + + finish: + + free(cumnr); + + free(ix_aux); + + free(tmin); + + free(tmax); + + free(dt); + + free(r2_all); + + return exitflag; + + failed: + + if (exitflag == -1) + + f_ap_print_error (__FUNCTION__, line, AP_ERR_MSG_MEM); + + goto finish; + +} + + + + +int f_ap_s_Ub_init (const double* s, const double* Ub, long* ix, const int D, \ + const long* N, const long ns, double* out_s, double* out_Ub) +{ +/* S H O R T D E S C R I P T I O N + * + * Remaps elements of the signal, modulator upper bound, and interpolation index + * arrays to comply with the Intel's MKL DFT indexing convention. */ + +/* I N P U T A R G U M E N T S + * + * [s] - input signal. + * + * [Ub] - upper bound on the modulator. This array must have the same number of + * elements as the input signal. + * + * [ix] - array with the mapping between the indexes of the original and (possibly) + * interpolated signal arrays. This array contains ns number of elements. + * + * [D] - number of dimensions of the signal array. + * + * [N] - numbers of elements of the signal array in every dimension. + * + * [ns] - number of elements in the array ix. + */ + +/* O U T P U T A R G U M E N T S + * + * [exitflag] - exit flag. + * + * [out_s#] - output signal with rearranged placement and +2 elements in the last + * dimension to meet the MKL DFT requirements. + * + * [out_Ub#] - output upper bound on the modulator with rearranged placement and +2 + * elements in the last dimension to meet the MKL DFT requirements. + * + * [ix#] - array with the mapping between the indexes of the original and (possibly) + * interpolated signal arrays compliant with the Intel's MKL DFT indexing + * convention. + */ + +/* N O N S T A N D A R D F U N C T I O N S U S E D + * + * (1) f_ap_print_error. + */ + + + /* Declarations and initializations */ + + int line; + + int exitflag = 0; + + + long i0, i1, i2; + + long n = 1, n_2; + + long i_lin0, i_lin; + + long strd0[3], strd[3]; + + long *ilin0_to_ilin = NULL; + + + + /* Number of elements in the s and Ub arrays */ + + for (i0=0; i0 1) + { + ilin0_to_ilin = (long*) malloc(n*sizeof(long)); + + if (ilin0_to_ilin == NULL) {line=__LINE__-2; exitflag=-1; goto failed;} + } + + + + /* Remapping of the s and Ub arrays */ + + if (D == 1) + { + for (i0=0; i0 1) + + for (i0=0; i00; i--) + { + cs[i] = cs[i+1] * N[i]; + + rs[i] = rs[i+1] * N[i]; + } + + cs[0] = 0; + + rs[0] = 0; + + + + /* Setting strides of the DFT descriptor for the forward transform */ + + status = DftiSetValue(*dft_handle, DFTI_INPUT_STRIDES, rs); + + if (status != DFTI_NO_ERROR) {line=__LINE__-2; exitflag=-2; goto failed;} + + + status = DftiSetValue(*dft_handle, DFTI_OUTPUT_STRIDES, cs); + + if (status != DFTI_NO_ERROR) {line=__LINE__-2; exitflag=-2; goto failed;} + + + + /* Setting the calculations to be done in-place */ + + status = DftiSetValue (*dft_handle, DFTI_PLACEMENT, DFTI_INPLACE); + + if (status != DFTI_NO_ERROR) {line=__LINE__-2; exitflag=-2; goto failed;} + + + + /* Setting the data storage scheme in the Fourier domain */ + + status = DftiSetValue (*dft_handle, DFTI_CONJUGATE_EVEN_STORAGE, \ + DFTI_COMPLEX_COMPLEX); + + if (status != DFTI_NO_ERROR) {line=__LINE__-3; exitflag=-2; goto failed;} + + + status = DftiSetValue (*dft_handle, DFTI_PACKED_FORMAT, DFTI_CCE_FORMAT); + + if (status != DFTI_NO_ERROR) {line=__LINE__-2; exitflag=-2; goto failed;} + + + + /* Setting the backward transform as the inverse transform */ + + status = DftiSetValue (*dft_handle, DFTI_BACKWARD_SCALE, 1.0/((double)n)); + + if (status != DFTI_NO_ERROR) {line=__LINE__-2; exitflag=-2; goto failed;} + + + + /* Setting all computations to use one CPU thread */ + + status = DftiSetValue (*dft_handle, DFTI_THREAD_LIMIT, 1); + + if (status != DFTI_NO_ERROR) {line=__LINE__-2; exitflag=-2; goto failed;} + + + + /* Commitof the DFT descriptor */ + + status = DftiCommitDescriptor (*dft_handle); + + if (status != DFTI_NO_ERROR) {line=__LINE__-2; exitflag=-2; goto failed;} + + + + /* Output & Memory deallocation */ + + goto finish; + + finish: + + free(N_); + + free(rs); + + free(cs); + + return exitflag; + + failed: + + DftiFreeDescriptor(dft_handle); + + *dft_handle = 0; + + if (exitflag == -1) + + f_ap_print_error (__FUNCTION__, line, AP_ERR_MSG_MEM); + + else if (exitflag == -2) + + f_ap_print_error (__FUNCTION__, line, DftiErrorMessage(status)); + + goto finish; + +} + + + + +int f_ap_mkl_dft_PMw (double* s, const int D, const long* N, const long* iL, \ + const long* iR, DFTI_DESCRIPTOR_HANDLE* dft_handle) +{ +/* S H O R T D E S C R I P T I O N + * + * Implements the projection onto the set Mw by using the Intel's Mkl DFT routine. */ + +/* I N P U T A R G U M E N T S + * + * [s] - input signal + 2 additional array elements along the last dimension. + * + * [D] - number of dimensions of the signal array. + * + * [N] - numbers of elements of the signal array in every dimension. The additional + * two elements in the last dimension of s are not counted here. + * + * [iL] - indexes of the left cutoff frequencies. + * + * [iR] - indexes of the right cutoff frequencies. + * + * [dft_handle] - address of an empty variable for the comitted descriptor handle. + */ + +/* O U T P U T A R G U M E N T S + * + * [exitflag] - exit flag. + * + * [s#] - pointer to the updated and comitted descriptor handle. + */ + +/* N O N S T A N D A R D F U N C T I O N S U S E D + * + * (1) f_ap_print_error, (2) DftiSetValue, (3) DftiCommitDescriptor, + * + * (4) DftiComputeForward, (5) DftiComputeBackward. + */ + + + /* Declarations and initializations */ + + int line = 0; + + int exitflag = 0; + + + long i1, i2, i3; + + MKL_LONG status; + + MKL_LONG *rs = NULL; + + MKL_LONG *cs = NULL; + + + + /* Calculation */ + + if (D == 1) + { + /* Forward FFT */ + + status = DftiComputeForward (*dft_handle, s); + + if (status != DFTI_NO_ERROR) {line=__LINE__-2; exitflag=-2; goto failed;} + + + + /* Projection onto Mw in the Fourier domain */ + + for (i1 = 2*iL[0]; i1 < N[0]+2; i1++) + + s[i1] = 0; + + + + /* Backward FFT */ + + status = DftiComputeBackward (*dft_handle, s); + + if (status != DFTI_NO_ERROR) {line=__LINE__-2; exitflag=-2; goto failed;} + } + + else + { + /* Calculation of strides in real and conjugate domains */ + + rs = (MKL_LONG*) malloc((D+1)*sizeof(MKL_LONG)); + + if (rs == NULL) {line=__LINE__-2; exitflag=-1; goto failed;} + + + cs = (MKL_LONG*) malloc((D+1)*sizeof(MKL_LONG)); + + if (cs == NULL) {line=__LINE__-2; exitflag=-1; goto failed;} + + + cs[D] = 1; + + rs[D] = 1; + + cs[D-1] = (N[D-1]/2+1); + + rs[D-1] = cs[D-1]*2; + + for (i1=D-2; i1>0; i1--) + { + cs[i1] = cs[i1+1] * N[i1]; + + rs[i1] = rs[i1+1] * N[i1]; + } + + cs[0] = 0; + + rs[0] = 0; + + + + /* Update of strides of the FFT descriptor for the forward transform */ + + status = DftiSetValue(*dft_handle, DFTI_INPUT_STRIDES, rs); + + if (status != DFTI_NO_ERROR) {line=__LINE__-2; exitflag=-2; goto failed;} + + + status = DftiSetValue(*dft_handle, DFTI_OUTPUT_STRIDES, cs); + + if (status != DFTI_NO_ERROR) {line=__LINE__-2; exitflag=-2; goto failed;} + + + + /* Commit of the descriptor */ + + status = DftiCommitDescriptor (*dft_handle); + + if (status != DFTI_NO_ERROR) {line=__LINE__-2; exitflag=-2; goto failed;} +/*printf("\n\n"); +for (i1 = 0; i1 < N[0]; i1++) +{ + for (i2 = 0; i2 < N[1]+2; i2++) + + printf("%e ",s[rs[1]*i1+rs[2]*i2]); + + printf("\n\n"); +}*/ +/*mexErrMsgIdAndTxt("AP:Demodulation","Here we leave!");*/ + + + /* Forward FFT */ + + status = DftiComputeForward (*dft_handle, s); + + if (status != DFTI_NO_ERROR) {line=__LINE__-2; exitflag=-2; goto failed;} +/*printf("\n\n"); +for (i1 = 0; i1 < N[0]; i1++) +{ + for (i2 = 0; i2 < N[1]+2; i2++) + + printf("%e ",s[rs[1]*i1+rs[2]*i2]); + + printf("\n\n"); +}*/ + +/*printf("iL=%ld,iR=%ld\n",iL[0],iR[0]);*/ + /* Projection onto Mw in the Fourier domain */ + + if (D == 2) + { + for (i2 = 2*iL[1]; i2 < N[1]+2-(N[0]%2); i2++) + + for (i1 = 0; i1 < N[0]; i1++) + + s[rs[1]*i1+rs[2]*i2] = 0; + + for (i2 = 0; i2 < 2*iL[1]; i2++) + + for (i1 = iL[0]; i1 <= iR[0]; i1++) + + s[rs[1]*i1+rs[2]*i2] = 0; + } + + else if (D == 3) + { + for (i3 = 2*iL[2]; i3 < N[2]+2-(N[0]%2); i3++) + + for (i2 = 0; i2 < N[1]; i2++) + + for (i1 = 0; i1 < N[0]; i1++) + + s[rs[1]*i1+rs[2]*i2+rs[3]*i3] = 0; + + for (i3 = 0; i3 < 2*iL[2]; i3++) + + for (i2 = iL[1]; i2 <= iR[1]; i2++) + + for (i1 = 0; i1 < N[0]; i1++) + + s[rs[1]*i1+rs[2]*i2+rs[3]*i3] = 0; + + for (i3 = 0; i3 < 2*iL[2]; i3++) + + for (i1 = iL[0]; i1 <= iR[0]; i1++) + { + for (i2 = 0; i2 < iL[1]; i2++) + + s[rs[1]*i1+rs[2]*i2+rs[3]*i3] = 0; + + for (i2 = iR[1]+1; i2 < N[1]; i2++) + + s[rs[1]*i1+rs[2]*i2+rs[3]*i3] = 0; + } + } + + else + { + exitflag = 1; + + goto failed; + } +/*printf("\n\n"); +for (i1 = 0; i1 < N[0]; i1++) +{ + for (i2 = 0; i2 < N[1]+2; i2++) + + printf("%e ",s[rs[1]*i1+rs[2]*i2]); + + printf("\n\n"); +}*/ +/*mexErrMsgIdAndTxt("AP:Demodulation","Here we leave!");*/ + + /* Update of strides of the FFT descriptor for the backward transform */ + + status = DftiSetValue(*dft_handle, DFTI_INPUT_STRIDES, cs); + + if (status != DFTI_NO_ERROR) {line=__LINE__-2; exitflag=-2; goto failed;} + + + status = DftiSetValue(*dft_handle, DFTI_OUTPUT_STRIDES, rs); + + if (status != DFTI_NO_ERROR) {line=__LINE__-2; exitflag=-2; goto failed;} + + + + /* Commit of the descriptor */ + + status = DftiCommitDescriptor (*dft_handle); + + if (status != DFTI_NO_ERROR) {line=__LINE__-2; exitflag=-2; goto failed;} + + + + /* Backward FFT */ + + status = DftiComputeBackward (*dft_handle, s); + + if (status != DFTI_NO_ERROR) {line=__LINE__-2; exitflag=-2; goto failed;} + /*printf("\n\n"); +for (i1 = 0; i1 < N[0]; i1++) +{ + for (i2 = 0; i2 < N[1]+2; i2++) + + printf("%e ",s[rs[1]*i1+rs[2]*i2]); + + printf("\n\n"); +}*/ +/*mexErrMsgIdAndTxt("AP:Demodulation","Here we leave!");*/ + } + + + + /* Output & Memory deallocation */ + + goto finish; + + finish: + + free(rs); + + free(cs); + + return exitflag; + + failed: + + if (exitflag == -1) + + f_ap_print_error (__FUNCTION__, line, AP_ERR_MSG_MEM); + + else if (exitflag == 1) + + f_ap_print_error (__FUNCTION__, line, AP_ERR_MSG_D); + + else if (exitflag == -2) + + f_ap_print_error (__FUNCTION__, line, DftiErrorMessage(status)); + + goto finish; + +} + diff --git a/C/library/l_ap_input_validation.c b/C/library/l_ap_input_validation.c new file mode 100644 index 0000000..5149684 --- /dev/null +++ b/C/library/l_ap_input_validation.c @@ -0,0 +1,333 @@ + +/* C O P Y R I G H T N O T I C E + * + * Copyright ©2021. Institute of Science and Technology Austria (IST Austria). + * All Rights Reserved. The underlying technology is protected by PCT Patent + * Application No. PCT/EP2021/054650. + * + * This file is part of the AP demodulation library, which is free software: you can + * redistribute it and/or modify it under the terms of the GNU General Public License + * as published by the Free Software Foundation in version 2. + * + * This program 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 v2 for more details. You + * should have received a copy of the GNU General Public License v2 along with this + * program. If not, see https://www.gnu.org/licenses/. + * + * Contact the Technology Transfer Office, IST Austria, Am Campus 1, + * A-3400 Klosterneuburg, Austria, +43-(0)2243 9000, twist@ist.ac.at, for commercial + * licensing opportunities. + * + * All other inquiries should be directed to the author, Mantas Gabrielaitis, + * mgabriel@ist.ac.at + */ + + +/* C O N T E N T S + * + * A function for checking the validity of input parameters for f_ap_demodulation: + * + * (1) l_ap_input_validation. + */ + + + +int f_ap_input_validation (const double* s, const struct s_ParamsAP* Par, \ + const double* Ub, const double* t) +{ +/* SHORT DESCRIPTION + * + * Checks the validity of input parameters for f_ap_demodulation. + */ + +/* INPUT ARGUMENTS + * + * [s] - input signal. + * + * [Par] - pointer to the structure with demodulation parameters: + * + * .Al - demodulation algorithm. Possible options are: 'B' -- AP-Basic, + * 'A' -- AP-Accelerated, 'P' -- AP-Projected. + * + * .D - number of signal dimensions. + * + * .Fs - sampling frequencies. This is an array of D elements with sampling + * frequency for each dimension of the signal. + * + * .Fc - cutoff frequencies of the modulator. This is an array of D elements + * with cutoff frequencies for each dimension of the signal. + * + * .Et - infeasibility error tolerance used to control the termination of the + * demodulation algorithm. The iterative process is stopped when the + * infeasibility error, ϵ, drops to the level of .Et or below. If + * .Et ≤ 0, then the maximum allowed number of iterations, .Ni (see + * below), is completed. + * + * .Ni - maximum number of allowed iterations of the chosen AP algorithm. The + * chosen AP algorithm is iterated not more than .Ni times independent + * of whether ϵ drops at or below .Et. + * + * .Ns - number of elements in every dimension of the original input signal + * (before any possible interpolation). This is an array with the + * number of elements equal to the number of signal dimensions. It does + * not assume the additional two elements in the last dimension of s. + * + * .Nr - number of sample points on the refined uniform grid used for + * interpolating nonuniformly sampled signals. This is an array with + * the number of elements equal to the number of signal dimensions. + * This field must be defined only if the third input argument t is + * provided (see below). If t is not provided, .Nr is ignored. + * + * .Cp - compression parameter. If .Cp > 1, demodulation is performed by + * using signal compression (see equations (12)-(13) in the paper cited + * in the header of this file. This field is optional (the default is + * .Cp=1). + * + * .Br - indicator of premature termination of the 'AP-Accelerated' algorithm + * when the value of the λ factor drops below one. If .Br=1, premature + * termination is assumed (default). Otherwise, the AP-A is not stopped + * premature even when λ decreases below 1. This field is required only + * if .Al='A'. It is optional (the default is .Br=1). + * + * .im - array with the iteration numbers at which the modulator estimates + * have to be saved for the output. the first element is the length of + * the array (excluding the first element itself). + * + * .ie - array with the iteration numbers at which the infeasibility error, + * ϵ, values have to be saved for the output. The first element of .ie + * is the length of the array (excluding the first element itself). + * + * [Ub] - upper bound on the modulator. This array must have the same number of + * elements as the input signal or must be set to NULL (if no upper bound on + * the modulator is assumed). + * + * [t] - sampling coordinates of the input signal. This is a 2D array with the number + * of columns equal to the dimension of the input signal and the number of rows + * equal to the number of signal sample points. This input argument is optional + * and should be used only if the input signal is sampled nonuniformly. + * Otherwise, t must be set to NULL. If t is not NULL, then the input argument + * Par must contain the field .Nr (see above). + */ + +/* O U T P U T A R G U M E N T S + * + * [exitflag] - exit flag (+1 indicates an error). + */ + + + + /* Declaration and initialization of the local variables */ + + int exitflag = 0; + + long i; + + long Ns = 1; + + long Ns3; + + + + + /* Input validation */ + + if (Par->Al != 'B' && Par->Al != 'A' && Par->Al != 'P') + { + + f_ap_print_error (__FUNCTION__, __LINE__, AP_ERR_MSG_AL); + + exitflag = 1; + } + + else if (Par->D<1 || Par->D>3) + { + f_ap_print_error (__FUNCTION__, __LINE__, AP_ERR_MSG_D); + + exitflag = 1; + } + + else if ( (Par->D>0 && (Par->Fs[0] <= 0 || !isfinite(Par->Fs[0]))) || \ + (Par->D>1 && (Par->Fs[1] <= 0 || !isfinite(Par->Fs[1]))) || \ + (Par->D>2 && (Par->Fs[2] <= 0 || !isfinite(Par->Fs[2]))) ) + { + f_ap_print_error (__FUNCTION__, __LINE__, AP_ERR_MSG_FS); + + exitflag = 1; + } + + else if ( (Par->D>0 && (Par->Fc[0] <= 0 || !isfinite(Par->Fc[0]))) || \ + (Par->D>1 && (Par->Fc[1] <= 0 || !isfinite(Par->Fc[1]))) || \ + (Par->D>2 && (Par->Fc[2] <= 0 || !isfinite(Par->Fc[2]))) ) + { + f_ap_print_error (__FUNCTION__, __LINE__, AP_ERR_MSG_FC); + + exitflag = 1; + } + + else if ( (Par->D>0 && Par->Fc[0] / Par->Fs[0] > 0.5) || \ + (Par->D>1 && Par->Fc[1] / Par->Fs[1] > 0.5) || \ + (Par->D>2 && Par->Fc[2] / Par->Fs[2] > 0.5) ) + { + f_ap_print_error (__FUNCTION__, __LINE__, AP_ERR_MSG_FC2); + + exitflag = 1; + } + + else if (!isfinite(Par->Et)) + { + f_ap_print_error (__FUNCTION__, __LINE__, AP_ERR_MSG_ET); + + exitflag = 1; + } + + else if (Par->Ni <= 0) + { + f_ap_print_error (__FUNCTION__, __LINE__, AP_ERR_MSG_NI); + + exitflag = 1; + } + + else if ( (t == NULL && \ + ((Par->D>0 && Par->Ns[0] <= 1) || \ + (Par->D>1 && Par->Ns[1] <= 1) || \ + (Par->D>2 && Par->Ns[2] <= 1))) || (t != NULL && Par->Ns[0] <= 1) ) + { + f_ap_print_error (__FUNCTION__, __LINE__, AP_ERR_MSG_NS); + + exitflag = 1; + } + + else if ( t != NULL && \ + ((Par->D>0 && Par->Nr[0] <= 1) || \ + (Par->D>1 && Par->Nr[1] <= 1) || \ + (Par->D>2 && Par->Nr[2] <= 1)) ) + { + f_ap_print_error (__FUNCTION__, __LINE__, AP_ERR_MSG_NR); + + exitflag = 1; + } + + else if (Par->Cp<1 || !isfinite(Par->Cp)) + { + f_ap_print_error (__FUNCTION__, __LINE__, AP_ERR_MSG_CP); + + exitflag = 1; + } + + else if (Par->Al == 'A' && Par->Br != 0 && Par->Br != 1) + { + f_ap_print_error (__FUNCTION__, __LINE__, AP_ERR_MSG_BR); + + exitflag = 1; + } + + else if (Par->ie[0] <= 0) + { + f_ap_print_error (__FUNCTION__, __LINE__, AP_ERR_MSG_IE); + + exitflag = 1; + } + + else if (Par->im[0] <= 0) + { + f_ap_print_error (__FUNCTION__, __LINE__, AP_ERR_MSG_IM); + + exitflag = 1; + } + + else + { + for (i = 1; i <= Par->ie[0]; i++) + { + if (Par->ie[i] < 0 || (i>1 && Par->ie[i-1] >= Par->ie[i])) + { + f_ap_print_error (__FUNCTION__, __LINE__, AP_ERR_MSG_IE2); + + exitflag = 1; + + break; + } + } + + + for (i = 1; i <= Par->im[0]; i++) + { + if (Par->im[i] < 0 || (i > 1 && Par->im[i-1] >= Par->im[i])) + { + f_ap_print_error (__FUNCTION__, __LINE__, AP_ERR_MSG_IM2); + + exitflag = 1; + + break; + } + } + } + + + + if (t == NULL) + { + for (i=0; i<(Par->D); i++) + + Ns = Ns * Par->Ns[i]; /* original signal */ + } + + else + + Ns = Par->Ns[0]; + + + + if (exitflag == 0) + { + for (i=0; i < Ns; i++) + { + if ( !isfinite(s[i]) ) + { + f_ap_print_error (__FUNCTION__, __LINE__, AP_ERR_MSG_S); + + exitflag = 1; + } + } + } + + + + if (exitflag == 0 && Ub != NULL) + { + + for (i=0; i < Ns; i++) + { + if ( !isfinite(Ub[i]) || Ub[i] < fabs(s[i])) + { + f_ap_print_error (__FUNCTION__, __LINE__, AP_ERR_MSG_UB); + + exitflag = 1; + } + } + } + + + + if (exitflag == 0 && t != NULL) + { + + Ns3 = Ns * Par->D; + + for (i=0; i < Ns3; i++) + { + if ( !isfinite(t[i]) ) + { + f_ap_print_error (__FUNCTION__, __LINE__, AP_ERR_MSG_T); + + exitflag = 1; + } + } + } + + + + return exitflag; + +} \ No newline at end of file diff --git a/COPYRIGHT b/COPYRIGHT new file mode 100644 index 0000000..026d536 --- /dev/null +++ b/COPYRIGHT @@ -0,0 +1,20 @@ +Copyright ©2021. Institute of Science and Technology Austria (IST Austria). +All Rights Reserved. The underlying technology is protected by PCT Patent Application +No. PCT/EP2021/054650. + +This file is part of the AP demodulation library, which is free software: you can +redistribute it and/or modify it under the terms of the GNU General Public License +as published by the Free Software Foundation in version 2. + +This program 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 v2 for more details. You +should have received a copy of the GNU General Public License v2 along with this +program. If not, see https://www.gnu.org/licenses/. + +Contact the Technology Transfer Office, IST Austria, Am Campus 1, +A-3400 Klosterneuburg, Austria, +43-(0)2243 9000, twist@ist.ac.at, for commercial +licensing opportunities. + +All other inquiries should be directed to the author, Mantas Gabrielaitis, +mgabriel@ist.ac.at \ No newline at end of file diff --git a/LICENSE b/LICENSE new file mode 100644 index 0000000..d159169 --- /dev/null +++ b/LICENSE @@ -0,0 +1,339 @@ + GNU GENERAL PUBLIC LICENSE + Version 2, June 1991 + + Copyright (C) 1989, 1991 Free Software Foundation, Inc., + 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA + Everyone is permitted to copy and distribute verbatim copies + of this license document, but changing it is not allowed. + + Preamble + + The licenses for most software are designed to take away your +freedom to share and change it. By contrast, the GNU General Public +License is intended to guarantee your freedom to share and change free +software--to make sure the software is free for all its users. This +General Public License applies to most of the Free Software +Foundation's software and to any other program whose authors commit to +using it. (Some other Free Software Foundation software is covered by +the GNU Lesser General Public License instead.) You can apply it to +your programs, too. + + When we speak of free software, we are referring to freedom, not +price. Our General Public Licenses are designed to make sure that you +have the freedom to distribute copies of free software (and charge for +this service if you wish), that you receive source code or can get it +if you want it, that you can change the software or use pieces of it +in new free programs; and that you know you can do these things. + + To protect your rights, we need to make restrictions that forbid +anyone to deny you these rights or to ask you to surrender the rights. +These restrictions translate to certain responsibilities for you if you +distribute copies of the software, or if you modify it. + + For example, if you distribute copies of such a program, whether +gratis or for a fee, you must give the recipients all the rights that +you have. You must make sure that they, too, receive or can get the +source code. And you must show them these terms so they know their +rights. + + We protect your rights with two steps: (1) copyright the software, and +(2) offer you this license which gives you legal permission to copy, +distribute and/or modify the software. + + Also, for each author's protection and ours, we want to make certain +that everyone understands that there is no warranty for this free +software. If the software is modified by someone else and passed on, we +want its recipients to know that what they have is not the original, so +that any problems introduced by others will not reflect on the original +authors' reputations. + + Finally, any free program is threatened constantly by software +patents. We wish to avoid the danger that redistributors of a free +program will individually obtain patent licenses, in effect making the +program proprietary. To prevent this, we have made it clear that any +patent must be licensed for everyone's free use or not licensed at all. + + The precise terms and conditions for copying, distribution and +modification follow. + + GNU GENERAL PUBLIC LICENSE + TERMS AND CONDITIONS FOR COPYING, DISTRIBUTION AND MODIFICATION + + 0. This License applies to any program or other work which contains +a notice placed by the copyright holder saying it may be distributed +under the terms of this General Public License. The "Program", below, +refers to any such program or work, and a "work based on the Program" +means either the Program or any derivative work under copyright law: +that is to say, a work containing the Program or a portion of it, +either verbatim or with modifications and/or translated into another +language. (Hereinafter, translation is included without limitation in +the term "modification".) Each licensee is addressed as "you". + +Activities other than copying, distribution and modification are not +covered by this License; they are outside its scope. The act of +running the Program is not restricted, and the output from the Program +is covered only if its contents constitute a work based on the +Program (independent of having been made by running the Program). +Whether that is true depends on what the Program does. + + 1. You may copy and distribute verbatim copies of the Program's +source code as you receive it, in any medium, provided that you +conspicuously and appropriately publish on each copy an appropriate +copyright notice and disclaimer of warranty; keep intact all the +notices that refer to this License and to the absence of any warranty; +and give any other recipients of the Program a copy of this License +along with the Program. + +You may charge a fee for the physical act of transferring a copy, and +you may at your option offer warranty protection in exchange for a fee. + + 2. You may modify your copy or copies of the Program or any portion +of it, thus forming a work based on the Program, and copy and +distribute such modifications or work under the terms of Section 1 +above, provided that you also meet all of these conditions: + + a) You must cause the modified files to carry prominent notices + stating that you changed the files and the date of any change. + + b) You must cause any work that you distribute or publish, that in + whole or in part contains or is derived from the Program or any + part thereof, to be licensed as a whole at no charge to all third + parties under the terms of this License. + + c) If the modified program normally reads commands interactively + when run, you must cause it, when started running for such + interactive use in the most ordinary way, to print or display an + announcement including an appropriate copyright notice and a + notice that there is no warranty (or else, saying that you provide + a warranty) and that users may redistribute the program under + these conditions, and telling the user how to view a copy of this + License. (Exception: if the Program itself is interactive but + does not normally print such an announcement, your work based on + the Program is not required to print an announcement.) + +These requirements apply to the modified work as a whole. If +identifiable sections of that work are not derived from the Program, +and can be reasonably considered independent and separate works in +themselves, then this License, and its terms, do not apply to those +sections when you distribute them as separate works. But when you +distribute the same sections as part of a whole which is a work based +on the Program, the distribution of the whole must be on the terms of +this License, whose permissions for other licensees extend to the +entire whole, and thus to each and every part regardless of who wrote it. + +Thus, it is not the intent of this section to claim rights or contest +your rights to work written entirely by you; rather, the intent is to +exercise the right to control the distribution of derivative or +collective works based on the Program. + +In addition, mere aggregation of another work not based on the Program +with the Program (or with a work based on the Program) on a volume of +a storage or distribution medium does not bring the other work under +the scope of this License. + + 3. You may copy and distribute the Program (or a work based on it, +under Section 2) in object code or executable form under the terms of +Sections 1 and 2 above provided that you also do one of the following: + + a) Accompany it with the complete corresponding machine-readable + source code, which must be distributed under the terms of Sections + 1 and 2 above on a medium customarily used for software interchange; or, + + b) Accompany it with a written offer, valid for at least three + years, to give any third party, for a charge no more than your + cost of physically performing source distribution, a complete + machine-readable copy of the corresponding source code, to be + distributed under the terms of Sections 1 and 2 above on a medium + customarily used for software interchange; or, + + c) Accompany it with the information you received as to the offer + to distribute corresponding source code. (This alternative is + allowed only for noncommercial distribution and only if you + received the program in object code or executable form with such + an offer, in accord with Subsection b above.) + +The source code for a work means the preferred form of the work for +making modifications to it. For an executable work, complete source +code means all the source code for all modules it contains, plus any +associated interface definition files, plus the scripts used to +control compilation and installation of the executable. However, as a +special exception, the source code distributed need not include +anything that is normally distributed (in either source or binary +form) with the major components (compiler, kernel, and so on) of the +operating system on which the executable runs, unless that component +itself accompanies the executable. + +If distribution of executable or object code is made by offering +access to copy from a designated place, then offering equivalent +access to copy the source code from the same place counts as +distribution of the source code, even though third parties are not +compelled to copy the source along with the object code. + + 4. You may not copy, modify, sublicense, or distribute the Program +except as expressly provided under this License. Any attempt +otherwise to copy, modify, sublicense or distribute the Program is +void, and will automatically terminate your rights under this License. +However, parties who have received copies, or rights, from you under +this License will not have their licenses terminated so long as such +parties remain in full compliance. + + 5. You are not required to accept this License, since you have not +signed it. However, nothing else grants you permission to modify or +distribute the Program or its derivative works. These actions are +prohibited by law if you do not accept this License. Therefore, by +modifying or distributing the Program (or any work based on the +Program), you indicate your acceptance of this License to do so, and +all its terms and conditions for copying, distributing or modifying +the Program or works based on it. + + 6. Each time you redistribute the Program (or any work based on the +Program), the recipient automatically receives a license from the +original licensor to copy, distribute or modify the Program subject to +these terms and conditions. You may not impose any further +restrictions on the recipients' exercise of the rights granted herein. +You are not responsible for enforcing compliance by third parties to +this License. + + 7. If, as a consequence of a court judgment or allegation of patent +infringement or for any other reason (not limited to patent issues), +conditions are imposed on you (whether by court order, agreement or +otherwise) that contradict the conditions of this License, they do not +excuse you from the conditions of this License. If you cannot +distribute so as to satisfy simultaneously your obligations under this +License and any other pertinent obligations, then as a consequence you +may not distribute the Program at all. For example, if a patent +license would not permit royalty-free redistribution of the Program by +all those who receive copies directly or indirectly through you, then +the only way you could satisfy both it and this License would be to +refrain entirely from distribution of the Program. + +If any portion of this section is held invalid or unenforceable under +any particular circumstance, the balance of the section is intended to +apply and the section as a whole is intended to apply in other +circumstances. + +It is not the purpose of this section to induce you to infringe any +patents or other property right claims or to contest validity of any +such claims; this section has the sole purpose of protecting the +integrity of the free software distribution system, which is +implemented by public license practices. Many people have made +generous contributions to the wide range of software distributed +through that system in reliance on consistent application of that +system; it is up to the author/donor to decide if he or she is willing +to distribute software through any other system and a licensee cannot +impose that choice. + +This section is intended to make thoroughly clear what is believed to +be a consequence of the rest of this License. + + 8. If the distribution and/or use of the Program is restricted in +certain countries either by patents or by copyrighted interfaces, the +original copyright holder who places the Program under this License +may add an explicit geographical distribution limitation excluding +those countries, so that distribution is permitted only in or among +countries not thus excluded. In such case, this License incorporates +the limitation as if written in the body of this License. + + 9. The Free Software Foundation may publish revised and/or new versions +of the General Public License from time to time. Such new versions will +be similar in spirit to the present version, but may differ in detail to +address new problems or concerns. + +Each version is given a distinguishing version number. If the Program +specifies a version number of this License which applies to it and "any +later version", you have the option of following the terms and conditions +either of that version or of any later version published by the Free +Software Foundation. If the Program does not specify a version number of +this License, you may choose any version ever published by the Free Software +Foundation. + + 10. If you wish to incorporate parts of the Program into other free +programs whose distribution conditions are different, write to the author +to ask for permission. For software which is copyrighted by the Free +Software Foundation, write to the Free Software Foundation; we sometimes +make exceptions for this. Our decision will be guided by the two goals +of preserving the free status of all derivatives of our free software and +of promoting the sharing and reuse of software generally. + + NO WARRANTY + + 11. BECAUSE THE PROGRAM IS LICENSED FREE OF CHARGE, THERE IS NO WARRANTY +FOR THE PROGRAM, TO THE EXTENT PERMITTED BY APPLICABLE LAW. EXCEPT WHEN +OTHERWISE STATED IN WRITING THE COPYRIGHT HOLDERS AND/OR OTHER PARTIES +PROVIDE THE PROGRAM "AS IS" WITHOUT WARRANTY OF ANY KIND, EITHER EXPRESSED +OR IMPLIED, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF +MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. THE ENTIRE RISK AS +TO THE QUALITY AND PERFORMANCE OF THE PROGRAM IS WITH YOU. SHOULD THE +PROGRAM PROVE DEFECTIVE, YOU ASSUME THE COST OF ALL NECESSARY SERVICING, +REPAIR OR CORRECTION. + + 12. IN NO EVENT UNLESS REQUIRED BY APPLICABLE LAW OR AGREED TO IN WRITING +WILL ANY COPYRIGHT HOLDER, OR ANY OTHER PARTY WHO MAY MODIFY AND/OR +REDISTRIBUTE THE PROGRAM AS PERMITTED ABOVE, BE LIABLE TO YOU FOR DAMAGES, +INCLUDING ANY GENERAL, SPECIAL, INCIDENTAL OR CONSEQUENTIAL DAMAGES ARISING +OUT OF THE USE OR INABILITY TO USE THE PROGRAM (INCLUDING BUT NOT LIMITED +TO LOSS OF DATA OR DATA BEING RENDERED INACCURATE OR LOSSES SUSTAINED BY +YOU OR THIRD PARTIES OR A FAILURE OF THE PROGRAM TO OPERATE WITH ANY OTHER +PROGRAMS), EVEN IF SUCH HOLDER OR OTHER PARTY HAS BEEN ADVISED OF THE +POSSIBILITY OF SUCH DAMAGES. + + END OF TERMS AND CONDITIONS + + How to Apply These Terms to Your New Programs + + If you develop a new program, and you want it to be of the greatest +possible use to the public, the best way to achieve this is to make it +free software which everyone can redistribute and change under these terms. + + To do so, attach the following notices to the program. It is safest +to attach them to the start of each source file to most effectively +convey the exclusion of warranty; and each file should have at least +the "copyright" line and a pointer to where the full notice is found. + + + Copyright (C) + + This program is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + + This program 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 General Public License along + with this program; if not, write to the Free Software Foundation, Inc., + 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. + +Also add information on how to contact you by electronic and paper mail. + +If the program is interactive, make it output a short notice like this +when it starts in an interactive mode: + + Gnomovision version 69, Copyright (C) year name of author + Gnomovision comes with ABSOLUTELY NO WARRANTY; for details type `show w'. + This is free software, and you are welcome to redistribute it + under certain conditions; type `show c' for details. + +The hypothetical commands `show w' and `show c' should show the appropriate +parts of the General Public License. Of course, the commands you use may +be called something other than `show w' and `show c'; they could even be +mouse-clicks or menu items--whatever suits your program. + +You should also get your employer (if you work as a programmer) or your +school, if any, to sign a "copyright disclaimer" for the program, if +necessary. Here is a sample; alter the names: + + Yoyodyne, Inc., hereby disclaims all copyright interest in the program + `Gnomovision' (which makes passes at compilers) written by James Hacker. + + , 1 April 1989 + Ty Coon, President of Vice + +This General Public License does not permit incorporating your program into +proprietary programs. If your program is a subroutine library, you may +consider it more useful to permit linking proprietary applications with the +library. If this is what you want to do, use the GNU Lesser General +Public License instead of this License. diff --git a/MATLAB/README.txt b/MATLAB/README.txt new file mode 100644 index 0000000..a428e29 --- /dev/null +++ b/MATLAB/README.txt @@ -0,0 +1,90 @@ + + +C O N T E N T S O F T H E C U R R E N T D I R E C T O R Y + +[./library_m] - folder with the MATLAB m-file code of the AP demodulation library: + + (i) f_ap_demodulation.m - defines the f_ap_demodulation function (the user's + interface to the AP demodulation algorithms). + + (ii) f_ap_basic.m - defines the f_ap_basic function, which estimates the + modulator of a signal by using the AP-Basic algorithm. + + (iii) f_ap_accelerated.m - defines the f_ap_accelerated function, which + estimates the modulator of a signal by using the + AP-Accelerated algorithm. + + (iv) f_ap_projected.m - defines the f_ap_accelerated function, which estimates + the modulator of a signal by using the AP-Projected + algorithm. + + +[./library_mex] - folder with the code of the MATLAB MEX function interface of the + AP demodulation library implemented in C: + + (i) f_ap_demodulation_mex.c - A MEX gateway function for f_ap_demodulation + defined in ../C/library/. + + (ii) compilation.m - MATLAB script for compiling f_ap_demodulation_mex.c to + produce the corresponding MEX function executable. + + +[./examples] - folder with five examples (example[1-5].m) of signal demodulation, + demonstrating various usage cases of f_ap_demodulation and + f_ap_demodulation_mex. + + + +U S A G E + +The AP demodulation approach was formulated and developed in + + M. Gabrielaitis. "Fast and Accurate Amplitude Demodulation of Wideband Signals," + IEEE Transactions on Signal Processing, vol. 69, pp. 4039-4054, 2021. DOI: + 10.1109/TSP.2021.3087899. + +We refer the user to this publication for any questions related to the working +principle of the new method. + +All demodulation capabilities of the current library are accessible solely via the +f_ap_demodulation (MATLAB m-file function) or f_ap_demodulation_mex (MATLAB MEX +function executable). The sufficient information for the correct use of these +functions can be found in the comments section of, respectively, f_ap_demodulation.m +and f_ap_demodulation_mex.c. Specifically, the requirements for and properties of the +input and output arguments are described there. To consolidate the confidence in +using f_ap_demodulation, consider analyzing the five examples of its applications +provided in ./examples. + + + +E X T E R N A L L I B R A R I E S (M E X F U N C T I O N) + +Besides the standard C library, generating the f_ap_demodulation_mex executable +requires Intel's oneAPI Math Kernel Library (oneMKL), specifically, its Fast Fourier +Transform routines. Thus, the latter must be installed on your system. oneMKL is +available for all major OS types free of charge (see https://software.intel.com/ +content/www/us/en/develop/tools/oneapi/base-toolkit/download.html for the +installation files). + + + +C O M P I L A T I O N A N D L I N K I N G F O R M E X F U N C T I O N + +f_ap_demodulation_mex.c is supposed to be compiled with GCC. This is the state-of- +the-art compiler available free of charge for all major OS types. It is preinstalled +on Linux systems as a rule. A GCC installation guide for Windows and Mac can be found +following this link: https://www.guru99.com/c-gcc-install.html. f_ap_demodulation_mex +executbale is generated with the MATLAB's built-in function mex. The latter must be +configured to use GCC for C language compilation (see MATLAB help for the +instructions). + +The compilation of f_ap_demodulation_mex.c is performed by running the compilation.m +script. The user, however, must set modify paths corresponding to the locations of +relevant AP demodulation and MKL library files on the system explicitly in the +mentioned script. This script is well commented on, making the required modification +straightforward. + +Note that the MKL library is dynamically linked. That is, the f_ap_demodulation_mex +executable cannot run successfully without the relevant MKL library files present on +the system during runtime. + diff --git a/MATLAB/examples/example1.m b/MATLAB/examples/example1.m new file mode 100644 index 0000000..7c6de53 --- /dev/null +++ b/MATLAB/examples/example1.m @@ -0,0 +1,193 @@ + +% C O P Y R I G H T N O T I C E +% +% Copyright ©2021. Institute of Science and Technology Austria (IST Austria). +% All Rights Reserved. The underlying technology is protected by PCT Patent +% Application No. PCT/EP2021/054650. +% +% This file is part of the AP demodulation library, which is free software: you can +% redistribute it and/or modify it under the terms of the GNU General Public License +% as published by the Free Software Foundation in version 2. +% +% This program 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 v2 for more details. You +% should have received a copy of the GNU General Public License v2 along with this +% program. If not, see https://www.gnu.org/licenses/. +% +% Contact the Technology Transfer Office, IST Austria, Am Campus 1, +% A-3400 Klosterneuburg, Austria, +43-(0)2243 9000, twist@ist.ac.at, for commercial +% licensing opportunities. +% +% All other inquiries should be directed to the author, Mantas Gabrielaitis, +% mgabriel@ist.ac.at + + + +% EXAMPLE 1 +% +% In this example, a synthetic 1D amplitude-modulated signal built from a harmonic +% carrier and a sinusoidal modulator is generated and demodulated by using the +% AP-Basic algorithm. The inferred modulator and carriers are then compared with the +% predefined ones. This example illustrates application of the functions +% 'f_ap_demodulation' and 'f_ap_demodulation_mex' in its simplest format, when the +% signal is sampled uniformly, and when only the final estimates of the modulator and +% infeasibility error are required. + + +close all + +clear all + + + +% Function type (m or mex) selected by the user + +fType = 'm'; % select between 'm' and 'mex' + + + +% Paths to the m and mex functions that perform AP-demodulation + +addpath ../library_m ../library_mex + + + +% Number of sample points + +n = 2^14; + + + +% Time vector + +T = linspace(0,25,n)'; + + + +% Modulator (nonnegative sinusoidal) + +m = (1.01+cos(2*pi*T)) / 2.01; + + + +% Carrier (harmonic) + +c = zeros(n,1); + +w = [0.4170, 0.7203, 0.0001, 0.3023, 0.1468, 0.0923, 0.1863, 0.3456, 0.3968, ... + 0.5388, 0.4192, 0.6852, 0.2045, 0.8781, 0.0274, 0.6705, 0.4173, 0.5587, ... + 0.1404, 0.1981, 0.8007, 0.9683, 0.3134, 0.6923, 0.8764, 0.8946, 0.0850, ... + 0.0391, 0.1698, 0.8781, 0.0983, 0.4211, 0.9579, 0.5332, 0.6919, 0.3155, ... + 0.6865, 0.8346, 0.0183, 0.7501]; + +for i=1:10 + + c = c + w(1+4*(i-1)) * cos(2*pi* (256*i*(0:n-1)'/n + w(2+4*(i-1)))); + + c = c + 0.01*w(3+4*(i-1)) * cos(2*pi* ((256*i+128)*(0:n-1)'/n + w(4*i))); + +end + +c = c / max(abs(c)); + + + +% Signal + +s = m .* c; + + + +% Demodulation + +if strcmpi(fType, 'm') + + fDemod = @(x,y) f_ap_demodulation (x, y); + +elseif strcmpi(fType, 'mex') + + fDemod = @(x,y) f_ap_demodulation_mex (x, y); + +end + + +Par.Al = 'B'; + +Par.Fs = 1/(T(2)-T(1)); + +Par.Fc = 1.5; + +Par.Et = 10^-5; + +Par.Ni = 10^3; + + +[m_, e, iter, tcpu] = fDemod (s, Par); + + + +% Visualization: signal and modulator + +figure('Name', 'Signal and Modulator') + +set(gcf, 'Units', 'centimeters'); + +set(gcf, 'Position', [8 17 35 7]); + + +plot(T,abs(s),'LineWidth',0.75) + +hold on + +plot(T,m,'k','LineWidth',2) + +plot(T,m_,'r','LineWidth',1) + +hold off + + +xlabel('t') + +ylabel('x') + +axis([3.5 8.5 0 1]) + +box('off') + + +legend({'$\bf{|s|}$','$\bf{m}$','$\hat{\bf{m}}$'}, 'Interpreter','latex', ... + 'Location',[0.545, 0.75, 0.1, 0.1], 'FontSize',11) + + + +% Visualization: carrier + +figure('Name', 'Carrier') + +set(gcf, 'Units', 'centimeters'); + +set(gcf, 'Position', [8 6 35 7]); + + +plot(T,c,'k','LineWidth',2) + +hold on + +plot(T,s./m_,'Color',[0 0.7 0],'LineWidth',0.75) + +hold off + + +xlabel('t') + +ylabel('c') + +axis([3.5 4.0 -1 1]) + +box('off') + + +legend({'$\bf{c}$','$\hat{\bf{c}}$'}, 'Interpreter','latex', ... + 'Location',[0.55, 0.80, 0.1, 0.1], 'FontSize',11) + diff --git a/MATLAB/examples/example2.m b/MATLAB/examples/example2.m new file mode 100644 index 0000000..a4f3ab9 --- /dev/null +++ b/MATLAB/examples/example2.m @@ -0,0 +1,245 @@ + +% C O P Y R I G H T N O T I C E +% +% Copyright ©2021. Institute of Science and Technology Austria (IST Austria). +% All Rights Reserved. The underlying technology is protected by PCT Patent +% Application No. PCT/EP2021/054650. +% +% This file is part of the AP demodulation library, which is free software: you can +% redistribute it and/or modify it under the terms of the GNU General Public License +% as published by the Free Software Foundation in version 2. +% +% This program 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 v2 for more details. You +% should have received a copy of the GNU General Public License v2 along with this +% program. If not, see https://www.gnu.org/licenses/. +% +% Contact the Technology Transfer Office, IST Austria, Am Campus 1, +% A-3400 Klosterneuburg, Austria, +43-(0)2243 9000, twist@ist.ac.at, for commercial +% licensing opportunities. +% +% All other inquiries should be directed to the author, Mantas Gabrielaitis, +% mgabriel@ist.ac.at + + + +% EXAMPLE 2 +% +% In this example, a synthetic 2D amplitude-modulated signal built of a random-peaks +% carrier and an LP-random-field modulator is generated and demodulated by using the +% AP-Accelerated algorithm. The inferred modulator is then compared with the +% predefined one. This example illustrates an application of the functions +% 'f_ap_demodulation' and 'f_ap_demodulation_mex' for demodulating 2D signals. + + +close all + +clear all + + + +% Function type (m or mex) selected by the user + +fType = 'mex'; % select between 'm' and 'mex' + + + +% Paths to the m and mex functions that perform AP-demodulation + +addpath ../library_m ../library_mex + + + +% Sampling grid + +[X, Y] = meshgrid(0:0.005:1); % grid used for demodulation purposes + +n = size(X,1); + +[X2, Y2] = meshgrid(0:0.025:1); % grid used for visualization purposes + + + +% Modulator (LP-random-field) + +m = zeros(size(X)); + +w = [0.5173, 0.9470, 0.7655, 0.2824, 0.2210, 0.6862, 0.1671, 0.3924, 0.6181, ... + 0.4119, 0.0025, 0.8840]; + +j = 0; + +for i1=0:3 + + for i2=0:2 + + j = j + 1; + + m = m + w(j)*cos(2*pi*(i1*(0:n-1)/n+i2*(0:n-1)'/n)+w(j)); + + end + +end + +m = m - min(m(:)); + +m = m / max(m(:)); + + + +% Subsampled modulator + +m2 = interp2(X,Y,m,X2,Y2); + + + +% Carrier (random peaks) + +cc = [0.2220, 0.2067, 0.4884, 0.7659, 0.2968, 0.0807, 0.4413, 0.8799, 0.4142, ... +0.6288, 0.5999, 0.2847, 0.3276, 0.1656, 0.9602, 0.0243, 0.6998, 0.0229, ... +0.0016, 0.6398, 0.2591, 0.8705, 0.0022, 0.9815, 0.8137, 0.0291, 0.1115, ... +0.9649, 0.6354, 0.9267, 0.8248, 0.3610, 0.5464, 0.3655, 0.7951, 0.6389, ... +0.5835, 0.9435, 0.8436, 0.1008, 0.5104, 0.4783, 0.5147, 0.8005, 0.5726, ... +0.9851, 0.4524, 0.3320, 0.4077, 0.3303, 0.5267, 0.8930, 0.7699, 0.7100, ... +0.7673, 0.2396, 0.3636, 0.0601, 0.8132, 0.0634, 0.0851, 0.1703, 0.8146, ... +0.0598, 0.5710, 0.8257, 0.5809, 0.6523, 0.4593, 0.4967, 0.6040, 0.2185, ... +0.3653, 0.5421, 0.1344, 0.6792, 0.9972, 0.7268, 0.5885, 0.2456, 0.1656, ... +0.2782, 0.6248, 0.8938, 0.7703, 0.1198, 0.1547, 0.0004, 0.1611, 0.3383, ... +0.6337, 0.7140, 0.1907, 0.9386, 0.2925, 0.5712, 0.5882, 0.6689, 0.2065, ... +0.8112, 0.4473, 0.0477, 0.7143, 0.9132, 0.5831, 0.2310, 0.5443, 0.0688, ... +0.2748, 0.4697, 0.9382, 0.7194, 0.6532, 0.4846, 0.2203, 0.7646, 0.3682, ... +0.2952, 0.3210, 0.8742, 0.3407, 0.3807, 0.8874, 0.0466, 0.4245, 0.6895, ... +0.8714, 0.1656, 0.2115, 0.3212, 0.0040, 0.5278, 0.2544, 0.0021, 0.7132, ... +0.2993, 0.1577, 0.0797, 0.4320, 0.7166, 0.6996, 0.2169, 0.4249, 0.3518, ... +0.9538, 0.0420, 0.2370, 0.2778, 0.1176, 0.9295, 0.9049, 0.9100, 0.0045, ... +0.4728, 0.2697, 0.8605, 0.3467, 0.1765, 0.1566, 0.3896, 0.3706, 0.0925, ... +0.2739, 0.5697, 0.1154, 0.5252, 0.1085, 0.6850, 0.5321, 0.2710, 0.4101, ... +0.4250, 0.6209, 0.4344, 0.0603, 0.6348, 0.8840, 0.1613, 0.2017, 0.9767, ... +0.4458, 0.0140, 0.5626, 0.0921, 0.3696, 0.1786, 0.5415, 0.3791, 0.4891, ... +0.9199, 0.7476, 0.0325, 0.8198, 0.6506, 0.4079, 0.5986, 0.1102, 0.9745, ... +0.2820, 0.5881, 0.6075, 0.0562, 0.1923, 0.9325, 0.5379, 0.8779, 0.8283, ... +0.9809, 0.7651, 0.7090, 0.9867, 0.8185, 0.5057, 0.9999, 0.2444, 0.0221, ... +0.4203, 0.9045, 0.7554, 0.7669, 0.2311, 0.5881, 0.3348, 0.6998, 0.6015, ... +0.8316, 0.9592, 0.6635, 0.7464, 0.5475, 0.4010, 0.8478, 0.7038, 0.6468, ... +0.2885, 0.3513, 0.8803, 0.4530, 0.8347, 0.8748, 0.1201, 0.4327, 0.1652, ... +0.9487, 0.1985, 0.1338, 0.5014, 0.0863, 0.9117, 0.2444, 0.0224, 0.7622, ... +0.1201, 0.5540, 0.4990, 0.0593, 0.4552, 0.8290, 0.7764, 0.3111, 0.0139, ... +0.4933, 0.9800, 0.9281, 0.5271, 0.6275, 0.8195, 0.2567, 0.6348, 0.9209; ... +... +0.8707, 0.9186, 0.6117, 0.5184, 0.1877, 0.7384, 0.1583, 0.2741, 0.2961, ... +0.5798, 0.2658, 0.2536, 0.1442, 0.9639, 0.1884, 0.2046, 0.7795, 0.5777, ... +0.5155, 0.9856, 0.8025, 0.9227, 0.4695, 0.3989, 0.5465, 0.0865, 0.2512, ... +0.6318, 0.8119, 0.9126, 0.0942, 0.0355, 0.7961, 0.2443, 0.3521, 0.4934, ... +0.9393, 0.1117, 0.3460, 0.3834, 0.9611, 0.8500, 0.4466, 0.0204, 0.4114, ... +0.8014, 0.7029, 0.3600, 0.8986, 0.0827, 0.6608, 0.9652, 0.7591, 0.7016, ... +0.9743, 0.2215, 0.8103, 0.4497, 0.2642, 0.2421, 0.8078, 0.1953, 0.8103, ... +0.9650, 0.3025, 0.6594, 0.4728, 0.2419, 0.4187, 0.2111, 0.7486, 0.4256, ... +0.4779, 0.2652, 0.3019, 0.6013, 0.5609, 0.6158, 0.6050, 0.0713, 0.0604, ... +0.3483, 0.1597, 0.7306, 0.1815, 0.9303, 0.4071, 0.8455, 0.7393, 0.7428, ... +0.8889, 0.2651, 0.6342, 0.2582, 0.4256, 0.9901, 0.8879, 0.3714, 0.3307, ... +0.4458, 0.0275, 0.8636, 0.8304, 0.4334, 0.3488, 0.0260, 0.0195, 0.5254, ... +0.9706, 0.2852, 0.9667, 0.1369, 0.7670, 0.0974, 0.7489, 0.6573, 0.6528, ... +0.7640, 0.5043, 0.4854, 0.2949, 0.5953, 0.0584, 0.3418, 0.8476, 0.4459, ... +0.1661, 0.5777, 0.9752, 0.0141, 0.7735, 0.7273, 0.9084, 0.6823, 0.1977, ... +0.8807, 0.8074, 0.0984, 0.5922, 0.5683, 0.3182, 0.4868, 0.0836, 0.6991, ... +0.0567, 0.6351, 0.7000, 0.6762, 0.1547, 0.6705, 0.6147, 0.3651, 0.2509, ... +0.9189, 0.4966, 0.5290, 0.8573, 0.8691, 0.2597, 0.7464, 0.4243, 0.0416, ... +0.1331, 0.8380, 0.6698, 0.3229, 0.8539, 0.9468, 0.9107, 0.5435, 0.5347, ... +0.9536, 0.4167, 0.4788, 0.1567, 0.6897, 0.5744, 0.5319, 0.2763, 0.7413, ... +0.2140, 0.8967, 0.1410, 0.5868, 0.1654, 0.1175, 0.0989, 0.1180, 0.5080, ... +0.8384, 0.3963, 0.0145, 0.9733, 0.5382, 0.3626, 0.5241, 0.4665, 0.8944, ... +0.0441, 0.0947, 0.2165, 0.3988, 0.6922, 0.5553, 0.1900, 0.6840, 0.6099, ... +0.2820, 0.8589, 0.5188, 0.6836, 0.1999, 0.7759, 0.1299, 0.6388, 0.3019, ... +0.6396, 0.0041, 0.0733, 0.2318, 0.5796, 0.0011, 0.9319, 0.8927, 0.6560, ... +0.7422, 0.3567, 0.0286, 0.0101, 0.5152, 0.6869, 0.0101, 0.0947, 0.3179, ... +0.6055, 0.9869, 0.8656, 0.3513, 0.8832, 0.8192, 0.9975, 0.7987, 0.4752, ... +0.4678, 0.3848, 0.1052, 0.5652, 0.3117, 0.5121, 0.1686, 0.7248, 0.9140, ... +0.7688, 0.5681, 0.3609, 0.6887, 0.9962, 0.4963, 0.6052, 0.5642, 0.1326, ... +0.1580, 0.0161, 0.7922, 0.4014, 0.0568, 0.3983, 0.3021, 0.1127, 0.3214]'; + +c = zeros(size(X)); + +for i=1:size(cc,1) + + c = c + exp(-((X-cc(i,1)).^2+(Y-cc(i,2)).^2)*8000); + +end + + + +% Signal + +s = c .* m; + + + +% Demodulation + +if strcmpi(fType, 'm') + + fDemod = @(x,y) f_ap_demodulation (x, y); + +elseif strcmpi(fType, 'mex') + + fDemod = @(x,y) f_ap_demodulation_mex (x, y); + +end + + +Par.Al = 'A'; + +Par.Fs = [n n]; + +Par.Fc = [4 4]; + +Par.Et = 10^-6; + +Par.Ni = 10^2; + + +[m_, e, iter, tcpu] = fDemod (s, Par); + + + +% Visualization: signal and modulator + +figure('Name', 'Signal and Modulator') + +set(gcf, 'Units', 'centimeters'); + +set(gcf, 'Position', [15 17 20 15]); + + +surf(X,Y,s,'FaceColor','none', 'EdgeColor', [0 0 0], 'EdgeAlpha',0.4, ... + 'LineWidth',0.2) + +hold on + +colormap summer + +surf(X,Y,m_, 'EdgeColor','none', 'FaceColor', 'flat', 'FaceAlpha',0.6) + +surf(X2,Y2,m2, 'FaceColor','none', 'EdgeColor', [1 1 1]*0.85, 'EdgeAlpha',0.6, ... + 'LineWidth',1.25) + +hold off + + +xlabel('t_1') + +ylabel('t_2') + +zlabel('x') + + +view(-198,59) + +grid('on') + +box('on') + + +legend({'$\bf{|s|}$','$\bf{m}$','$\hat{\bf{m}}$'}, 'Interpreter','latex', ... + 'Location',[0.59, 0.80, 0.1, 0.1], 'FontSize',11) + diff --git a/MATLAB/examples/example3.m b/MATLAB/examples/example3.m new file mode 100644 index 0000000..5843da6 --- /dev/null +++ b/MATLAB/examples/example3.m @@ -0,0 +1,166 @@ + +% C O P Y R I G H T N O T I C E +% +% Copyright ©2021. Institute of Science and Technology Austria (IST Austria). +% All Rights Reserved. The underlying technology is protected by PCT Patent +% Application No. PCT/EP2021/054650. +% +% This file is part of the AP demodulation library, which is free software: you can +% redistribute it and/or modify it under the terms of the GNU General Public License +% as published by the Free Software Foundation in version 2. +% +% This program 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 v2 for more details. You +% should have received a copy of the GNU General Public License v2 along with this +% program. If not, see https://www.gnu.org/licenses/. +% +% Contact the Technology Transfer Office, IST Austria, Am Campus 1, +% A-3400 Klosterneuburg, Austria, +43-(0)2243 9000, twist@ist.ac.at, for commercial +% licensing opportunities. +% +% All other inquiries should be directed to the author, Mantas Gabrielaitis, +% mgabriel@ist.ac.at + + + +% EXAMPLE 3 +% +% In this example, a nonuniformly sampled 1D amplitude-modulated signal built from a +% white-noise carrier and an LP-random modulator is generated and demodulated by +% using the AP-Projected algorithm. The inferred modulator is then compared with the +% predefined one. This example illustrates application of the functions +% 'f_ap_demodulation' and 'f_ap_demodulation_mex' for demodulating nonuniformly +% sampled signals. + + +close all + +clear all + + + +% Function type (m or mex) selected by the user + +fType = 'mex'; % select between 'm' and 'mex' + + + +% Paths to the m and mex functions that perform AP-demodulation + +addpath ../library_m ../library_mex + + + +% Number of sample points + +n = 2^10; + + + +% Time vector + +rng(1); + +t = [0; cumsum(0.5+2.0*rand(n,1))]; + +t = 10 * t(1:end-1) / t(end); + + + +% Modulator (LP-random) + +m = zeros(n,1); + +w = [1.5648, 0.5312, 0.1413, 0.7588, -0.8616, -0.3586, 0.9106, -0.1787, -0.0108, ... + -0.0989, -0.3559, -0.4015, 0.2917, -0.3458, -1.1990, 0.7651, -0.9884, -1.1668, ... + 0.6584, -1.3693, -0.4143, 0.6981, 0.9181, 2.9477, -0.3757, 0.8218, -0.5699, ... + -0.5495, -0.0108, 0.0854, 1.0189, 1.7904, 0.7577, -0.8434, -0.7300, -0.8675, ... + -0.6159, -0.3319, 1.4449, 0.6741]; + +for i=0:19 + + m = m + w(1+2*i)*cos(2*pi*i*t/10+w(2*(i+1))); + +end + +m = m - min(m(:)) + 0.001; + +m = m / max(m(:)); + + + +% Carrier (sinusoidal) + +c = sin(2*pi*5*t); + + + +% Signal + +s = c .* m; + + + +% Demodulation + +if strcmpi(fType, 'm') + + fDemod = @(x,y,z) f_ap_demodulation (x, y, [], z); + +elseif strcmpi(fType, 'mex') + + fDemod = @(x,y,z) f_ap_demodulation_mex (x, y, [], z); + +end + + +Par.Al = 'P'; + +Par.Fs = 1 / mean(diff(t)); + +Par.Fc = 25 * Par.Fs / (n*8); + +Par.Et = 10^-4; + +Par.Ni = 3*10^4; + +Par.Nr = n * 8; + + +[m_, e, Niter, tcpu] = fDemod (s, Par, t); + + + +% Visualization: signal and modulator + +figure('Name', 'Signal and Modulator') + +set(gcf, 'Units', 'centimeters'); + +set(gcf, 'Position', [8 17 35 7]); + + +plot(t,abs(s),'o-','LineWidth',0.75,'MarkerSize',2) + +hold on + +plot(t,m,'k','LineWidth',2) + +plot(t,m_,'r','LineWidth',1) + +hold off + + +xlabel('t') + +ylabel('x') + +axis([2 4 0 1]) + +box('off') + + +legend({'$\bf{|s|}$','$\bf{m}$','$\hat{\bf{m}}$'}, 'Interpreter','latex', ... + 'Location',[0.40, 0.75, 0.1, 0.1], 'FontSize',11) + diff --git a/MATLAB/examples/example4.m b/MATLAB/examples/example4.m new file mode 100644 index 0000000..9d8d2c8 --- /dev/null +++ b/MATLAB/examples/example4.m @@ -0,0 +1,253 @@ + +% C O P Y R I G H T N O T I C E +% +% Copyright ©2021. Institute of Science and Technology Austria (IST Austria). +% All Rights Reserved. The underlying technology is protected by PCT Patent +% Application No. PCT/EP2021/054650. +% +% This file is part of the AP demodulation library, which is free software: you can +% redistribute it and/or modify it under the terms of the GNU General Public License +% as published by the Free Software Foundation in version 2. +% +% This program 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 v2 for more details. You +% should have received a copy of the GNU General Public License v2 along with this +% program. If not, see https://www.gnu.org/licenses/. +% +% Contact the Technology Transfer Office, IST Austria, Am Campus 1, +% A-3400 Klosterneuburg, Austria, +43-(0)2243 9000, twist@ist.ac.at, for commercial +% licensing opportunities. +% +% All other inquiries should be directed to the author, Mantas Gabrielaitis, +% mgabriel@ist.ac.at + + + +% EXAMPLE 4 +% +% In this example, a synthetic 1D amplitude-modulated signal built of a random-spikes +% carrier and an LP-random modulator is generated and demodulated by using the +% AP-Accelerated algorithm. Demodulation is performed with and without (default) +% upper bound constraints. The inferred modulators are then compared with the +% predefined one. Moreover, the intermediate estimates of the modulators and +% infeasibility errors are saved and then visualized. This example illustrates how to +% set the upper bound constraints on the modulator estimates and obtain intermediate +% modulator estimates and infeasibility errors by using functions 'f_ap_demodulation' +% and 'f_ap_demodulation_mex'. The results obtained also illustrate that imposing an +% upper bound may reduce the rate of convergence of the AP algorithm in terms of +% infeasibility error. However, that has no practical consequences on the convergence +% in terms of the demodulation error. + + +close all + +clear all + + + +% Function type (m or mex) selected by the user + +fType = 'm'; % select between 'm' and 'mex' + + + +% Paths to the m and mex functions that perform AP-demodulation + +addpath ../library_m ../library_mex + + +% Number of sample points + +n = 2^10; + + +% Time vector + +T = linspace(0,10,n)'; + + + +% Modulator (LP-random) + +m = zeros(n,1); + +w = [1.5648, 0.5312, 0.1413, 0.7588, -0.8616, -0.3586, 0.9106, -0.1787, -0.0108, ... + -0.0989, -0.3559, -0.4015, 0.2917, -0.3458, -1.1990, 0.7651, -0.9884, -1.1668, ... + 0.6584, -1.3693]; + +for i=0:9 + + m = m + w(1+2*i)*cos(2*pi*i*(0:n-1)'/n+w(2*(i+1))); + +end + +m = m - min(m(:)) + 0.001; + +m = m / max(m(:)); + + + +% Carrier (random-spikes) + +ci = [1, 38, 46, 27, 35, 30, 29, 32, 36, 37, 41, 38, 45, 32, 50, 27, 45, 38, 42, ... + 30, 32, 48, 53, 35, 45, 50, 51]; + +ci = cumsum(ci); + +c = zeros(n,1); + +c(ci) = 1; + + + +% Signal + +s = c .* m; + + + +% Demodulation (without UB) + +if strcmpi(fType, 'm') + + fDemod = @(x,y,z) f_ap_demodulation (x, y, z); + +elseif strcmpi(fType, 'mex') + + fDemod = @(x,y,z) f_ap_demodulation_mex (x, y, z); + +end + + +Par.Al = 'A'; + +Par.Fs = 1 / (T(2)-T(1)); + +Par.Fc = 10 * Par.Fs / n; + +Par.Et = -1; + +Par.Ni = 10^3; + +Par.im = 1:10^3; + +Par.ie = 1:10^3; + + +[m_1, e1] = fDemod (s, Par, []); + + + +% Demodulation (with UB) + +Ub = [0.59*ones(120,1); 0.27*ones(40,1); 0.59*ones(140,1); 1*ones(100,1); ... + 0.71*ones(623,1); 0.311]; + + +[m_2, e2] = fDemod (s, Par, Ub); + + + +% Demodulation errors + +E1 = zeros(numel(Par.im),1); + +E2 = zeros(numel(Par.im),1); + +for i=1:numel(Par.im) + + E1(i) = norm(m_1(:,i)-m) / norm(m); + + E2(i) = norm(m_2(:,i)-m) / norm(m); + +end + + + +% Visualization: signal and modulator + +figure('Name', 'Signal and Modulator') + +set(gcf, 'Units', 'centimeters'); + +set(gcf, 'Position', [8 17 35 7]); + + +plot(T,abs(s),'o-','LineWidth',0.75,'MarkerSize',2) + +hold on + +plot(T,m, 'k','LineWidth',2) + +plot(T,m_1(:,end), 'r','LineWidth',1.5) + +plot(T,m_2(:,end), 'Color',[0.18, 0.74, 0.19],'LineWidth',1.5) + +plot(T,Ub, '--','Color',0.6*[1 1 1],'LineWidth',1) + +hold off + + +xlabel('t') + +ylabel('x') + +% axis([2 4 0 1]) + +box('off') + + +legend({'$\bf{|s|}$','$\bf{m}$','$\hat{\bf{m}}$ wo ub','$\hat{\bf{m}}$ w ub','ub'}, ... + 'Interpreter','latex', 'Location',[0.50, 0.75, 0.1, 0.1], 'FontSize',11) + + + +% Visualization: errors + +figure('Name', 'Infeasibility and Demodulation Errors') + +set(gcf, 'Units', 'centimeters'); + +set(gcf, 'Position', [8 6 35 7]); + + +subplot(1,2,1) + +loglog(e1, 'r','LineWidth',1.5) + +hold on + +loglog(e2, 'Color',[0.18, 0.74, 0.19],'LineWidth',1.5) + +hold off + + +xlabel('N_{iter}') + +ylabel('\epsilon') + + +legend({'$\epsilon$ wo ub','$\epsilon$ w ub'}, ... + 'Interpreter','latex', 'Location',[0.30, 0.55, 0.1, 0.1], 'FontSize',11) + + +subplot(1,2,2) + +semilogx(E1, 'r','LineWidth',1.5) + +hold on + +semilogx(E2, 'Color',[0.18, 0.74, 0.19],'LineWidth',1.5) + +hold off + + +xlabel('N_{iter}') + +ylabel('E') + + +legend({'$E_m$ wo ub','$E_m$ w ub'}, ... + 'Interpreter','latex', 'Location',[0.75, 0.75, 0.1, 0.1], 'FontSize',11) + diff --git a/MATLAB/examples/example5.m b/MATLAB/examples/example5.m new file mode 100644 index 0000000..5e7eb42 --- /dev/null +++ b/MATLAB/examples/example5.m @@ -0,0 +1,200 @@ + +% C O P Y R I G H T N O T I C E +% +% Copyright ©2021. Institute of Science and Technology Austria (IST Austria). +% All Rights Reserved. The underlying technology is protected by PCT Patent +% Application No. PCT/EP2021/054650. +% +% This file is part of the AP demodulation library, which is free software: you can +% redistribute it and/or modify it under the terms of the GNU General Public License +% as published by the Free Software Foundation in version 2. +% +% This program 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 v2 for more details. You +% should have received a copy of the GNU General Public License v2 along with this +% program. If not, see https://www.gnu.org/licenses/. +% +% Contact the Technology Transfer Office, IST Austria, Am Campus 1, +% A-3400 Klosterneuburg, Austria, +43-(0)2243 9000, twist@ist.ac.at, for commercial +% licensing opportunities. +% +% All other inquiries should be directed to the author, Mantas Gabrielaitis, +% mgabriel@ist.ac.at + + + +% EXAMPLE 5 +% +% In this example, a synthetic 1D amplitude-modulated signal built from a +% regular-spikes carrier and two LP-random modulators is generated and demodulated by +% using the AP-Basic algorithm. The two modulators are used to shape the lower and +% upper envelopes of the signal. The inferred modulators are compared with the +% predefined ones. This example illustrates how to infer the upper and lower +% envelopes of a signal by using functions 'f_ap_demodulation' and +% 'f_ap_demodulation_mex'. + + +close all + +clear all + + + +% Function type (m or mex) selected by the user + +fType = 'm'; % select between 'm' and 'mex' + + + +% Paths to the m and mex functions that perform AP-demodulation + +addpath ../library_m ../library_mex + + + +% Number of sample points + +n = 2^10; + + + +% Time vector + +T = linspace(0,10,n)'; + + + +% Modulators (LP-random) + +m1 = zeros(n,1); + +w = [1.5648, 0.5312, 0.1413, 0.7588, -0.8616, -0.3586, 0.9106, -0.1787, -0.0108, ... + -0.0989, -0.3559, -0.4015, 0.2917, -0.3458, -1.1990, 0.7651, -0.9884, -1.1668, ... + 0.6584, -1.3693, 0.7608, 0.7810, 0.9041, 0.2338, 0.1767, 0.3911, 0.3206, ... + 0.8155, 0.6135, 0.7600]; + +for i=0:9 + + m1 = m1 + w(1+2*i)*cos(2*pi*i*(0:n-1)'/n+w(2*(i+1))); + +end + +m1 = m1 - min(m1(:)); + +m1 = m1 / max(m1(:)); + + +m2 = zeros(n,1); + +for i=14:-1:0 + + m2 = m2 + w(1+2*i)*cos(2*pi*i*(0:n-1)'/n+w(2*(i+1))); + +end + +m2 = m2(end:-1:1) - min(m2(:)); + +m2 = m2 / max(m2(:)) / 2; + + + +% Carrier (regular spikes) + +ci = 4:32:n; + +c1 = zeros(n,1); + +c2 = zeros(n,1); + +c1(ci) = 1; + +c2(ci-2) = -1; + +c2(ci+2) = -1; + +c = c1 + c2; + + + +% Signal + +s = c1 .* m1 + c2 .* m2; + + + +% Demodulation + +if strcmpi(fType, 'm') + + fDemod = @(x,y) f_ap_demodulation (x, y); + +elseif strcmpi(fType, 'mex') + + fDemod = @(x,y) f_ap_demodulation_mex (x, y); + +end + + +Par.Al = 'B'; + +Par.Fs = 1 / (T(2)-T(1)); + +Par.Fc = 15 * Par.Fs / n; + +Par.Et = 10^-6; + +Par.Ni = 10^3; + + +smmin = min(-s); + +[m1_, e1] = fDemod (-s-smmin, Par); + +m1_ = -m1_ - smmin; + +smin = min(s); + +[m2_, e2] = fDemod (s-smin, Par); + +m2_ = m2_ + smin; + + + +% Visualization: signal and modulator + +figure('Name', 'Signal and Modulators') + +set(gcf, 'Units', 'centimeters'); + +set(gcf, 'Position', [8 17 35 7]); + + +pl(1) = plot(T,s,'LineWidth',0.75,'MarkerSize',2); + +hold on + +pl(2) = plot(T,m1,'k','LineWidth',2); + +pl(3) = plot(T,-m2,'k','LineWidth',2); + +pl(4) = plot(T,m1_,'r--','LineWidth',2); + +pl(5) = plot(T,m2_,'--','Color',[0.18, 0.74, 0.19],'LineWidth',1); + +hold off + + +xlabel('t') + +ylabel('x') + +axis([0 10 -0.5 1]) + +box('off') + + +legend([pl(1:2) pl(4:5)], {'$\bf{|s|}$','$\bf{m}$ lower and upper',... + '$\hat{\bf{m}}$ lower','$\hat{\bf{m}}$ upper'}, 'Interpreter','latex', ... + 'Location',[0.53, 0.80, 0.1, 0.1], 'FontSize',11) + diff --git a/MATLAB/library_m/f_ap_accelerated.m b/MATLAB/library_m/f_ap_accelerated.m new file mode 100644 index 0000000..b3adf5a --- /dev/null +++ b/MATLAB/library_m/f_ap_accelerated.m @@ -0,0 +1,352 @@ + +% C O P Y R I G H T N O T I C E +% +% Copyright ©2021. Institute of Science and Technology Austria (IST Austria). +% All Rights Reserved. The underlying technology is protected by PCT Patent +% Application No. PCT/EP2021/054650. +% +% This file is part of the AP demodulation library, which is free software: you can +% redistribute it and/or modify it under the terms of the GNU General Public License +% as published by the Free Software Foundation in version 2. +% +% This program 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 v2 for more details. You +% should have received a copy of the GNU General Public License v2 along with this +% program. If not, see https://www.gnu.org/licenses/. +% +% Contact the Technology Transfer Office, IST Austria, Am Campus 1, +% A-3400 Klosterneuburg, Austria, +43-(0)2243 9000, twist@ist.ac.at, for commercial +% licensing opportunities. +% +% All other inquiries should be directed to the author, Mantas Gabrielaitis, +% mgabriel@ist.ac.at + + + +function [m_out, e_out, iter] = f_ap_accelerated (s, Par, Ub, ix) +% > +% S H O R T D E S C R I P T I O N +% +% Calculates the modulator of a signal by using the AP-Basic algorithm. Signals +% defined in up to 3 dimensions are allowed. +% +% +% I N P U T A R G U M E N T S +% +% [s] - input signal. This is a real array of a chosen dimension D ≤ 3. +% +% [Par] - parameters characterizing the signal and demodulation procedure. This is a +% variable of the structure type. Its fields are as follows: +% +% .Fs - sampling frequencies for each dimension of the signal. This is an +% array of D elements. +% +% .Fc - cutoff frequencies of the modulator for each dimension of the signal. +% This is an array of D elements. +% +% .Et - infeasibility error tolerance used to control the termination of +% the demodulation algorithm. The iterative process is stopped when the +% infeasibility error, ϵ, drops to the level of .Et or below. If +% .Et ≤ 0, then the maximum allowed number of iterations, .Ni (see +% below), is completed. +% +% .Ni - maximum number of allowed iterations of the chosen AP algorithm. The +% chosen AP algorithm is iterated not more than .Ni times independent +% of whether ϵ drops at or below .Et. +% +% .Ns numbers of sample points of the original input signal (before any +% possible interpolation) in every dimension. This is an array with the +% number of elements equal to the number of signal dimensions. +% +% .Nx - numbers of sample points of the provided input signal (possibly +% interpolated) in every dimension. This is an array with the number +% of elements equal to the number of signal dimensions. +% +% .Br - indicator of premature termination of the 'AP-Accelerated' algorithm +% when the value of the λ factor drops below one. If .Br=1, premature +% termination is assumed (default). Otherwise, the AP-A is not stopped +% premature even when λ decreases below 1. This field is required only +% if .Al='A'. It is optional (the default is .Br=1). +% +% .im - array with the iteration numbers at which the modulator estimates +% have to be saved for the output. If .im is empty, only the final +% modulator estimate is saved. The same is true when numel(i.im)=1 and +% .im(1)=.Ni. This field is optional (.im=[] is assumed by default). +% +% .ie - array with the iteration numbers at which the infeasibility error, ϵ, +% values have to be saved for the output. If .ie is empty, only the +% final error is saved. The same is true when numel(i.im)=1 and +% .im(1)=.Ni. This field is optional (.ie=[] is assumed by default). +% +% [Ub] - array with the upper bound constraints on the modulator. If Ub is an empty +% array, no upper bound constraint is assumed. Otherwise, Ub must be an array +% with the same dimensions as s. This input argument is optional (the default +% is Ub=[]). +% +% [ix] - linear indexes of the elements of an interpolated signal corresponding to +% every element of the original signal, i.e., mapping between the linear +% indexes of the original and interpolated signals. +% +% +% O U T P U T A R G U M E N T S +% +% [m_out] - modulator estimates at predefined interations (D+1 array). +% +% [e_out] - infeasibility error estimates at predefined iterations (1D array). +% +% [Ni_out] - number of iterations actually used (scalar). + + + +% [1] DEFINITION AND INITIALIZATION OF RELEVANT VARIABLES + +% Dimensionality of the input signal + +D = sum(size(s)~=1); + + + +% Numbers of sample points of the original and possibly interpolated signals + +ns = prod(Par.Ns); + +nx = prod(Par.Nx); + + + +% Full-wave rectification and scaling of the signal + +if D == 1 + + s = s(:); + +end + +s_abs = abs(s); % absolute value of the signal + +max_s = max(s_abs(:)); % maximum of the absolute-value signal + +s_abs = s_abs / max_s; % scaled absolute signal + +if ~isempty(Ub) + + Ub = Ub / max_s; + +end + + + +% Frequencies to filter out with the projection onto Mw + +fL = 1 + ceil(Par.Fc(:) ./ (Par.Fs(:)./Par.Nx)); + +fR = Par.Nx - fL + 2; + +i_cutoff = cell(D,1); + +for i=1:D + + i_cutoff{i} = fL(i)+1:fR(i)-1; + +end + + +% Initialization of the infeasibility error variable + +if Par.Et > 0 + + Etol = (Par.Et / max_s)^2; + +else + + Etol = Par.Et; + +end + +E = sum(s_abs(:).^2) / nx; + + + +% Initialization of the output arguments and related variables + +iter_m = 1; % iteration variable for storing the modulator estimates + +nim = numel(Par.im); + + +m_out = zeros([Par.Ns; max(numel(Par.im),1)]'); % output array with modulator estimates + +if ~isempty(Par.im) && Par.im(iter_m) == 1 + + m_out(1+(iter_m-1)*nx:iter_m*nx) = s_abs; + + iter_m = iter_m + 1; + +end + + + +iter_e = 1; % iteration variable for storing the infeasibility errors + +nie = numel(Par.ie); + + +e_out = zeros(numel(Par.ie),1); % output array with the infeasibility errors + +if ~isempty(Par.ie) && Par.ie(iter_e) == 1 + + e_out(iter_e) = E; + + iter_e = iter_e + 1; + +end + + + +% [3] CALCULATIONS + +iter = 0; + +m = s_abs; + +a = zeros(size(m)); + +nom = sum(m(:).*m(:)); + + +while E > Etol && iter < Par.Ni + + iter = iter + 1; + + + % Projection onto Mw + + b = fftn(m-a); + + if D == 1 + + b(i_cutoff{1}) = 0; + + elseif D == 2 + + b(i_cutoff{1},:) = 0; + + b(:,i_cutoff{2}) = 0; + + elseif D == 3 + + b(i_cutoff{1},:,:) = 0; + + b(:,i_cutoff{2},:) = 0; + + b(:,:,i_cutoff{3}) = 0; + + end + + b = ifftn(b); + + + % Lambda + + denom = sum(b(:).*b(:)); + + if denom ~= 0 + + la = nom / denom; + + else + + la = 1; + + end + + if la < 1 && Par.Br == 1 + + break; + + end + + + % Updated a + + a = a + la * b; + + + % Projection onto Cd + + m = a; + + m(m < s_abs) = s_abs(m < s_abs); + + if ~isempty(Ub) + + m(m > Ub) = Ub(m > Ub); + + end + + + % Infeasibility error + + nom = sum((m(:)-a(:)).^2); + + E = nom / nx; + + + % Output (modulator) + + if iter_m <= nim && (iter == Par.im(iter_m) || (E <= Etol && nim == 1 && ... + Par.im(1) == Par.Ni)) + + m_out(1+(iter_m-1)*ns:iter_m*ns) = m(abs(ix)); + + iter_m = iter_m + 1; + + end + + + % Output (infeasibility error) + + if iter_e <= nie && (iter == Par.ie(iter_e) || (E <= Etol && nie == 1 && ... + Par.ie(1) == Par.Ni)) + + e_out(iter_e) = E; + + iter_e = iter_e + 1; + + end + +end + + + +% [4] FINAL OUTPUT + +% If only the final estimate of the modulator is requested + +if isempty(Par.im) + + m_out(:) = m(abs(ix)); + +end + + +% Scaling back the modulator estimates + +m_out = max_s * m_out; + + +if isempty(Par.ie) + + + e_out = E; + +end + + +% Scaling back the error estimates + +e_out = max_s * sqrt(e_out); + + + diff --git a/MATLAB/library_m/f_ap_basic.m b/MATLAB/library_m/f_ap_basic.m new file mode 100644 index 0000000..d89e636 --- /dev/null +++ b/MATLAB/library_m/f_ap_basic.m @@ -0,0 +1,314 @@ + +% C O P Y R I G H T N O T I C E +% +% Copyright ©2021. Institute of Science and Technology Austria (IST Austria). +% All Rights Reserved. The underlying technology is protected by PCT Patent +% Application No. PCT/EP2021/054650. +% +% This file is part of the AP demodulation library, which is free software: you can +% redistribute it and/or modify it under the terms of the GNU General Public License +% as published by the Free Software Foundation in version 2. +% +% This program 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 v2 for more details. You +% should have received a copy of the GNU General Public License v2 along with this +% program. If not, see https://www.gnu.org/licenses/. +% +% Contact the Technology Transfer Office, IST Austria, Am Campus 1, +% A-3400 Klosterneuburg, Austria, +43-(0)2243 9000, twist@ist.ac.at, for commercial +% licensing opportunities. +% +% All other inquiries should be directed to the author, Mantas Gabrielaitis, +% mgabriel@ist.ac.at + + + +function [m_out, e_out, iter] = f_ap_basic (s, Par, Ub, ix) +% > +% S H O R T D E S C R I P T I O N +% +% Calculates the modulator of a signal by using the AP-Basic algorithm. Signals +% defined in up to 3 dimensions are allowed. +% +% +% I N P U T A R G U M E N T S +% +% [s] - input signal. This is a real array of a chosen dimension D ≤ 3. +% +% [Par] - parameters characterizing the signal and demodulation procedure. This is a +% variable of the structure type. Its fields are as follows: +% +% .Fs - sampling frequencies for each dimension of the signal. This is an +% array of D elements. +% +% .Fc - cutoff frequencies of the modulator for each dimension of the signal. +% This is an array of D elements. +% +% .Et - infeasibility error tolerance used to control the termination of +% the demodulation algorithm. The iterative process is stopped when the +% infeasibility error, ϵ, drops to the level of .Et or below. If +% .Et ≤ 0, then the maximum allowed number of iterations, .Ni (see +% below), is completed. +% +% .Ni - maximum number of allowed iterations of the chosen AP algorithm. The +% chosen AP algorithm is iterated not more than .Ni times independent +% of whether ϵ drops at or below .Et. +% +% .Ns numbers of sample points of the original input signal (before any +% possible interpolation) in every dimension. This is an array with the +% number of elements equal to the number of signal dimensions. +% +% .Nx - numbers of sample points of the provided input signal (possibly +% interpolated) in every dimension. This is an array with the number +% of elements equal to the number of signal dimensions. +% +% .im - array with the iteration numbers at which the modulator estimates +% have to be saved for the output. If .im is empty, only the final +% modulator estimate is saved. The same is true when numel(i.im)=1 and +% .im(1)=.Ni. This field is optional (.im=[] is assumed by default). +% +% .ie - array with the iteration numbers at which the infeasibility error, ϵ, +% values have to be saved for the output. If .ie is empty, only the +% final error is saved. The same is true when numel(i.im)=1 and +% .im(1)=.Ni. This field is optional (.ie=[] is assumed by default). +% +% [Ub] - array with the upper bound constraints on the modulator. If Ub is an empty +% array, no upper bound constraint is assumed. Otherwise, Ub must be an array +% with the same dimensions as s. This input argument is optional (the default +% is Ub=[]). +% +% [ix] - linear indexes of the elements of an interpolated signal corresponding to +% every element of the original signal, i.e., mapping between the linear +% indexes of the original and interpolated signals. +% +% +% O U T P U T A R G U M E N T S +% +% [m_out] - modulator estimates at predefined interations (D+1 array). +% +% [e_out] - infeasibility error estimates at predefined iterations (1D array). +% +% [Ni_out] - number of iterations actually used (scalar). + + + +% [1] DEFINITION AND INITIALIZATION OF RELEVANT VARIABLES + +% Dimensionality of the input signal + +D = sum(size(s)~=1); + + + +% Numbers of sample points of the original and possibly interpolated signals + +ns = prod(Par.Ns); + +nx = prod(Par.Nx); + + + +% Full-wave rectification and scaling of the signal + +if D == 1 + + s = s(:); + +end + +s_abs = abs(s); % absolute value of the signal + +max_s = max(s_abs(:)); % maximum of the absolute-value signal + +s_abs = s_abs / max_s; % scaled absolute signal + +if ~isempty(Ub) + + Ub = Ub / max_s; + +end + + + +% Frequencies to filter out with the projection onto Mw + +fL = 1 + ceil(Par.Fc(:) ./ (Par.Fs(:)./Par.Nx)); + +fR = Par.Nx - fL + 2; + +i_cutoff = cell(D,1); + +for i=1:D + + i_cutoff{i} = fL(i)+1:fR(i)-1; + +end + + +% Initialization of the infeasibility error variable + +if Par.Et > 0 + + Etol = (Par.Et / max_s)^2; + +else + + Etol = Par.Et; + +end + +E = sum(s_abs(:).^2) / nx; + + + +% Initialization of the output arguments and related variables + +iter_m = 1; % iteration variable for storing the modulator estimates + +nim = numel(Par.im); + + +m_out = zeros([Par.Ns; max(numel(Par.im),1)]'); % output array with modulator estimates + +if ~isempty(Par.im) && Par.im(iter_m) == 1 + + m_out(1+(iter_m-1)*nx:iter_m*nx) = s_abs; + + iter_m = iter_m + 1; + +end + + + +iter_e = 1; % iteration variable for storing the infeasibility errors + +nie = numel(Par.ie); + + +e_out = zeros(numel(Par.ie),1); % output array with the infeasibility errors + +if ~isempty(Par.ie) && Par.ie(iter_e) == 1 + + e_out(iter_e) = E; + + iter_e = iter_e + 1; + +end + + + +% [3] CALCULATIONS + +iter = 0; + +m = s_abs; + + +while E > Etol && iter < Par.Ni + + iter = iter + 1; + + + % Projection onto Mw + + a = m; + + a = fftn(a); + + if D == 1 + + a(i_cutoff{1}) = 0; + + elseif D == 2 + + a(i_cutoff{1},:) = 0; + + a(:,i_cutoff{2}) = 0; + + elseif D == 3 + + a(i_cutoff{1},:,:) = 0; + + a(:,i_cutoff{2},:) = 0; + + a(:,:,i_cutoff{3}) = 0; + + end + + a = ifftn(a); + + + % Projection onto Cd + + m = a; + + m(m < s_abs) = s_abs(m < s_abs); + + if ~isempty(Ub) + + m(m > Ub) = Ub(m > Ub); + + end + + + % Infeasibility error + + E = sum((a(:)-m(:)).^2) / nx; + + + % Output (modulator) + + if iter_m <= nim && (iter == Par.im(iter_m) || (E <= Etol && nim == 1 && ... + Par.im(1) == Par.Ni)) + + m_out(1+(iter_m-1)*ns:iter_m*ns) = m(abs(ix)); + + iter_m = iter_m + 1; + + end + + + % Output (infeasibility error) + + if iter_e <= nie && (iter == Par.ie(iter_e) || (E <= Etol && nie == 1 && ... + Par.ie(1) == Par.Ni)) + + e_out(iter_e) = E; + + iter_e = iter_e + 1; + + end + +end + + + +% [4] FINAL OUTPUT + +% If only the final estimate of the modulator is requested + +if isempty(Par.im) + + m_out(:) = m(abs(ix)); + +end + + +% Scaling back the modulator estimates + +m_out = max_s * m_out; + + +if isempty(Par.ie) + + + e_out = E; + +end + + +% Scaling back the error estimates + +e_out = max_s * sqrt(e_out); + diff --git a/MATLAB/library_m/f_ap_demodulation.m b/MATLAB/library_m/f_ap_demodulation.m new file mode 100644 index 0000000..6f59d75 --- /dev/null +++ b/MATLAB/library_m/f_ap_demodulation.m @@ -0,0 +1,700 @@ + +% C O P Y R I G H T N O T I C E +% +% Copyright ©2021. Institute of Science and Technology Austria (IST Austria). +% All Rights Reserved. The underlying technology is protected by PCT Patent +% Application No. PCT/EP2021/054650. +% +% This file is part of the AP demodulation library, which is free software: you can +% redistribute it and/or modify it under the terms of the GNU General Public License +% as published by the Free Software Foundation in version 2. +% +% This program 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 v2 for more details. You +% should have received a copy of the GNU General Public License v2 along with this +% program. If not, see https://www.gnu.org/licenses/. +% +% Contact the Technology Transfer Office, IST Austria, Am Campus 1, +% A-3400 Klosterneuburg, Austria, +43-(0)2243 9000, twist@ist.ac.at, for commercial +% licensing opportunities. +% +% All other inquiries should be directed to the author, Mantas Gabrielaitis, +% mgabriel@ist.ac.at + + + +% C O N T E N T S +% +% f_ap_demodulation (main function), +% +% f_ap_input_validation (local function), +% +% f_ap_compression (local function), +% +% f_ap_interpolation (local function). + + + +function [m_out, e_out, iter, tcpu] = f_ap_demodulation (s, Par, Ub, t) +% > +% S H O R T D E S C R I P T I O N +% +% f_ap_demodulation performs demodulation of the input signal in the offline mode by +% using the AP demodulation approach formulated and developed in +% +% M. Gabrielaitis. "Fast and Accurate Amplitude Demodulation of Wideband signals," +% IEEE Transactions on Signal Processing, vol. 69, pp. 4039-4054, 2021. DOI: +% 10.1109/TSP.2021.3087899 +% +% Please cite the above publication when referring to the algorithms implemented in +% this code or the code itself. +% +% See https://github.com/mgabriel-lt/ap-demodulation for the latest version of the +% code and user-friendly explanations on the working principle, domains of +% application, and advice on the usage of different AP demodulation algorithms in +% practice. +% +% +% I N P U T A R G U M E N T S +% +% [s] - input signal. This is a real array of a chosen dimension D ≤ 3. If the input +% signal is nonuniformly sampled, and the coordinates of the sample points are +% provided via the fourth input argument, t (see below), s must be a 1D array, +% independent of the true dimensionality of the signal. +% +% [Par] - parameters characterizing the signal and demodulation procedure. This is a +% variable of the structure type. Its fields are as follows: +% +% .Al - demodulation algorithm. Possible options are: 'B' - AP-Basic, +% 'A' - AP-Accelerated, 'P' - AP-Projected. +% +% .Fs - sampling frequencies for each dimension of the signal. This is an +% array of D elements. +% +% .Fc - cutoff frequencies of the modulator for each dimension of the signal. +% This is an array of D elements. +% +% .Et - infeasibility error tolerance used to control the termination of +% the demodulation algorithm. The iterative process is stopped when the +% infeasibility error, ϵ, drops to the level of .Et or below. If +% .Et ≤ 0, then the maximum allowed number of iterations, .Ni (see +% below), is completed. +% +% .Ni - maximum number of allowed iterations of the chosen AP algorithm. The +% chosen AP algorithm is iterated not more than .Ni times independent +% of whether ϵ drops at or below .Et. +% +% .Nr - numbers of sample points on the refined uniform interpolation grid in +% every dimension. This is an array of D elements. This field must be +% defined only if the fourth input argument, t, is provided and is +% nonempy (see below). Otherwise, .Nr is ignored. +% +% .Cp - compression parameter. If .Cp > 1, demodulation is performed by using +% signal compression (see equations (12)-(13) in the paper cited in the +% header of this file. This field is optional (the default is .Cp=1). +% +% .Br - indicator of premature termination of the 'AP-Accelerated' algorithm +% when the value of the λ factor drops below one. If .Br=1, premature +% termination is assumed (default). Otherwise, the AP-A is not stopped +% premature even when λ decreases below 1. This field is required only +% if .Al='A'. It is optional (the default is .Br=1). +% +% .im - array with the iteration numbers at which the modulator estimates +% have to be saved for the output. If .im is empty, only the final +% modulator estimate is saved. The same is true when numel(i.im)=1 and +% .im(1)=.Ni. This field is optional (.im=[] is assumed by default). +% +% .ie - array with the iteration numbers at which the infeasibility error, ϵ, +% values have to be saved for the output. If .ie is empty, only the +% final error is saved. The same is true when numel(i.im)=1 and +% .im(1)=.Ni. This field is optional (.ie=[] is assumed by default). +% +% [Ub] - array with the upper bound constraints on the modulator. Two options are +% possible. If Ub is an empty array, no upper bound constraint is assumed. +% Otherwise, Ub must be an array with the same dimensions as s. This input +% argument is optional (Ub=[] is assumed by default). +% +% [t] - sampling coordinates. This is a real 2D array, with the number of rows equal +% to the number of signal sample points, and the number of columns equal to the +% dimension of the input signal. This input argument is optional and should +% be used only if the input signal is sampled nonuniformly. If t is provided, +% then the input argument Par must contain the field .Nr, and the input signal +% arrray, s, must be 1D, independent of the true dimensionality of the signal +% (see above). +% +% +% O U T P U T A R G U M E N T S +% +% [m_out] - modulator estimates at predefined interations (D+1 array). +% +% [e_out] - infeasibility error estimates at predefined iteratiHons (1D array). +% +% [iter] - number of iterations actually used (scalar). +% +% [tcpu] - running time of the function in seconds (scalar). + + +tic + +% [1] VALIDATION OF INPUT ARGUMENTS AND DEFAULT VALUES FOR UNDEFINED INPUTS + +if nargin < 3 + + Ub = []; + +end + +if nargin < 4 + + t = []; + +end + +if ~isfield(Par,'Cp') || isempty(Par.Cp) + + Par.Cp = 1; + +end + +if ~isfield(Par,'Br') || isempty(Par.Br) + + Par.Br = 1; + +end + +if ~isfield(Par,'im') + + Par.im = []; + +end + +if ~isfield(Par,'ie') + + Par.ie = []; + +end + +if nargin < 3 + + Ub = []; + +end + +f_ap_input_validation (s, Par, Ub, t) + + + +% [2] DEFINITION AND INITIALIZATION OF RELEVANT VARIABLES + +% Lengths of the original and possibly interpolated signal in each dimension + +if nargin > 3 && ~isempty(t) + + Par.Ns = numel(s); + + Par.Nx = Par.Nr; + +else + + Par.Ns = size(s)'; + + if sum(Par.Ns~=1) == 1 + + Par.Ns = max(Par.Ns); + + end + + Par.Nx = Par.Ns; + +end + + + +% [3] CALCULATIONS + +% Compression of the signal (if applicable) + +if Par.Cp > 1 + + s = f_ap_compression(s, 1/Par.Cp); + + if ~isempty(Ub) + + Ub = f_ap_compression(Ub, 1/Par.Cp); + + end + +end + + +% Interpolation of the signal (if applicable). + +if nargin > 3 && ~isempty(t) + + [s, Ub, ix] = f_ap_interpolation (s, Par, Ub, t); + +else + + ix = 1:numel(s); + +end + + +% Demodulation + +if Par.Al == 'B' + + [m_out, e_out, iter] = f_ap_basic (s, Par, Ub, ix); + +elseif Par.Al == 'A' + + [m_out, e_out, iter] = f_ap_accelerated (s, Par, Ub, ix); + +elseif Par.Al == 'P' + + [m_out, e_out, iter] = f_ap_projected (s, Par, Ub, ix); + +end + + +% Decompression + +if Par.Cp > 1 + + m_out = f_ap_compression(m_out, Par.Cp); % decompression of the signal + +end + + +% Computing time + +tcpu = toc; + +end + + + + +% LOCAL FUNCTIONS: f_ap_input_validation, f_ap_compression, f_ap_interpolation. + + +function f_ap_input_validation (s, Par, Ub, t) +% > +% S H O R T D E S C R I P T I O N +% +% Checks the validity of input arguments for f_ap_demodulation. +% +% +% I N P U T A R G U M E N T S +% +% [s] - input signal. This is a real array of a chosen dimension D ≤ 3. If the fourth +% input argument, t, is nonempty, i.e., interpolation of the input signal is +% required, then s must be a 1D array. +% +% [Par] - parameters characterizing the signal and demodulation procedure. This is a +% variable of the structure type. Its fields are as follows: +% +% .Al - demodulation algorithm. Possible options are: 'B' - AP-Basic, +% 'A' - AP-Accelerated, 'P' - AP-Projected. +% +% .Fs - sampling frequencies for each dimension of the signal. This is an +% array of D elements. +% +% .Fc - cutoff frequencies of the modulator for each dimension of the signal. +% This is an array of D elements. +% +% .Et - infeasibility error tolerance used to control the termination of +% the demodulation algorithm. The iterative process is stopped when the +% infeasibility error, ϵ, drops to the level of .Et or below. If +% .Et ≤ 0, then the maximum allowed number of iterations, .Ni (see +% below), is completed. +% +% .Ni - maximum number of allowed iterations of the chosen AP algorithm. The +% chosen AP algorithm is iterated not more than .Ni times independent +% of whether ϵ drops at or below .Et. +% +% .Nr - numbers of sample points on the refined uniform interpolation grid in +% every dimension. This is an array of D elements. This field must be +% defined only if the fourth input argument t is provided and is +% nonempy (see below). Otherwise, .Nr is ignored. +% +% .Cp - compression parameter. +% +% .Br - indicator of premature termination of the 'AP-Accelerated' algorithm +% when the value of the λ factor drops below one. +% +% .im - array with the iteration numbers at which the modulator estimates +% have to be saved for the output. If .im is empty, only the final +% modulator estimate is saved. +% +% .ie - array with the iteration numbers at which the infeasibility error, ϵ, +% values have to be saved for the output. If .ie is empty, only the +% final error is saved. +% +% [Ub] - array with the upper bound constraints on the modulator. This array must be +% empty or have the same dimensions as the first input argument, s. +% +% [t] - sampling coordinates. This is a real 2D array, with the number of rows equal +% to the number of signal sample points, and the number of columns equal to the +% dimension of the input signal. This input argument must be empty if no signal +% interpolation on a uniform grid is assumed. +% +% +% O U T P U T A R G U M E N T S +% +% None. + + + +% Dimensionality of the input signal + +if isempty(t) + + D = sum(size(s)>1); + +else + + D = size(t,2); + +end + + + +% Checking input argument s + +if numel(s) <= 1 || ndims(s) > 3 || ~isreal(s) + + error(['The array with the signal samples, s, must be real, not empty, and ',... + 'have no more than 3 dimensions!']) + +elseif ndims(s)==3 && any(size(s)==1) + + error(['If the array with the signal samples, s, is 3D, it must have no ',... + 'singleton dimensions!']) + +elseif ~isempty(t) && D > 1 + + error(['If the interpolation of the input signal, s, is assumed, s must be ',... + 'a 1D array!']) + +end + + + +% Checking fields of the input argument Par + +if ~all(isfield(Par, {'Fs','Fc','Al','Et','Ni'})) + + error(['The second input argument, Par, must be a structure with fields',... + '''Fs'', ''Fc'', ''Al'', ''Et'', and ''Ni'' at least!']) + +elseif ~isempty(t) && ~isfield(Par, 'Nr') + + error(['The second input argument, Par, must contain a field ''Nr'' with ',... + 'the size of the interpolated signal array if the fourth input argument ',... + 'is provided!']) + +end + + +% .Al +if Par.Al ~= 'B' && Par.Al ~= 'A' && Par.Al ~= 'P' && Par.Al ~= 'R' + + error(['The demodulation algorithm, Par.Al, must be one of ',... + '''B'', ''A'', ''P'', and ''R''!']); +% .Fs +elseif ~isreal(Par.Fs) || any(Par.Fs(:) <= 0) || numel(Par.Fs)~=D + + error(['The sampling frequency array , Par.Fs, must consist of positive ', ... + 'real entries whose number is equal to the number of dimensions of the ', ... + 'input signal array, s!']) +% .Fc +elseif ~isreal(Par.Fc) || any(Par.Fc(:) <= 0) || numel(Par.Fc)~=numel(Par.Fs) || ... + any(Par.Fc(:) > Par.Fs(:)*0.5) + + error(['The cutoff frequency array, Par.Fc, must consist of non-negative ', ... + 'real entries whose number is equal to the number of dimensions of the ', ... + 'input signal array, s. Moreover, Par.Fc must not exceed half of Par.Fs', ... + ' for each entry!']) +% .Et +elseif ~isreal(Par.Et) || ~isscalar(Par.Et) + + error('The infeasibility error tolerance, Par.Et, must be a real scalar!') +% .Ni +elseif Par.Ni < 0 || ~isreal(Par.Ni) || ~isscalar(Par.Ni) + + error(['The maximum number of allowed iterations, Par.Ni, must be a ', ... + 'non-negative real scalar!']) +% .Nr +elseif ~isempty(t) && (~isreal(Par.Nr) || any(Par.Nr <= 1) || ... + any(Par.Nr~=floor(Par.Nr)) || numel(Par.Nr)~=numel(Par.Fs)) + + error(['The array of numbers of sampling points, Par.Nr, must consist of ',... + 'integers exceeding 1 whose number is equal to the number of dimensions ',... + 'of the input signal array s!']) +% .Cp +elseif Par.Cp < 1 || ~isreal(Par.Cp) || ~isscalar(Par.Cp) + + error(['The compression parameter, Par.Cp, must be a real scalar not ', ... + 'smaller than 1!']) +% .Br +elseif ~isreal(Par.Br) || ~isscalar(Par.Br) + + error(['The indicator of the premature termination of the ''AP-Accelerated'','... + ' algorithm, Par.Br, must be a real scalar!']) +% .im +elseif any(Par.im <= 0) || ( ~isempty(Par.im) && ... + (~isreal(Par.im) || ~issorted(Par.im) || any(Par.im~=floor(Par.im)) || ... + numel(Par.im)~=numel(unique(Par.im))) ) + + error(['The array with the iteration numbers at which the modulator ', ... + 'estimates have to be saved for the output, Par.im, must be empty or ', ... + 'consist of unique nonnegative increasing integers!']) +% .ie +elseif any(Par.ie <= 0) || ( ~isempty(Par.ie) && ... + (~isreal(Par.ie) || ~issorted(Par.ie) || any(Par.ie~=floor(Par.ie)) || ... + numel(Par.ie)~=numel(unique(Par.ie))) ) + + error(['The array with the iteration numbers at which the modulator ', ... + 'estimates have to be saved for the output, Par.ie, must be empty or ', ... + 'consist of unique nonnegative increasing integers!']) + +end + + +% Checking input argument Ub + +if ~isempty(Ub) && ( ~isreal(Ub) || numel(Ub)~=numel(s) || any(Ub(:) +% S H O R T D E S C R I P T I O N +% +% Compresses the dynamic range of the input array by using the power function model. +% +% +% I N P U T A R G U M E N T S +% +% [in] - input signal. +% +% [Cp] - compression parameter. +% +% +% O U T P U T A R G U M E N T S +% +% [out] - compressed signal. + + + +out = sign(in) .* abs(in).^Cp; + + +end + + + + +function [out_s, out_Ub, out_ix] = f_ap_interpolation (s, Par, Ub, t) +% > +% S H O R T D E S C R I P T I O N +% +% Interpolates the input signal on a refined uniform grid following Eq. 23 in the +% paper cited in the header of the present file). +% +% +% I N P U T A R G U M E N T S +% +% [s] - input signal. +% +% [Par] - parameters characterizing the signal and demodulation procedure. This is a +% variable of the structure type. Its fields are as follows: +% +% .Nr - numbers of sample points on the refined uniform grid used for +% interpolating nonuniformly sampled signals. This is an array of D +% elements. +% +% [Ub] - array with the upper bound constraints on the modulator. This array must be +% empty or have the same number of elements as the input signal array, s. +% +% [t] - sampling coordinates. This is a real 2D array, with the number of rows equal +% to the number of signal sample points, and the number of columns equal to the +% dimension of the input signal. +% +% +% O U T P U T A R G U M E N T S +% +% [out_s] - interpolated signal. +% +% [out_Ub] - interpolated upper bounds on the modulator. +% +% [out_ix] - linear indexes of the elements of an interpolated signal corresponding +% to every element of the original signal, i.e., mapping between the +% linear indexes of the original and interpolated signals. + + + +% Overall length of the original signal + +ns = numel(s); + + + +% Cumulative product of the dimensions of the interpolated signal + +cumNr = cumprod(Par.Nr(1:end-1)); + +cumNr = cumNr(:); + + + +% New and old arrays of the signal and modulator's upper bound + +if numel(Par.Nr) == 1 + + out_s = zeros(Par.Nr,1); + +else + + out_s = zeros(Par.Nr); + +end + +if ~isempty(Ub) + + if numel(Par.Nr) == 1 + + out_Ub = zeros(Par.Nr,1) + Inf; + + else + + out_Ub = zeros(Par.Nr) + Inf; + + end + +else + + out_Ub = Ub; + +end + + + +% Bounds and step sizes of the time coordinates of the interpolated signal sample +% points + +t = t'; + +tmin = min(t,[],2); + +dt = (max(t,[],2) - tmin) ./ (Par.Nr(:) - 1); + + + +% Interpolation + +i_all = 1:ns; + +r_all = zeros(1,ns); + +out_ix = zeros(1,ns); + + +for i=1:ns + + % The closest point on the new grid + + ix = round((t(:,i)-tmin) ./ dt) + 1; + + + % The correponding distance + + r = norm((t(:,i)-tmin) - (ix-1).*dt); + + + % Linear index of the new point + + if sum(size(out_s)~=1) > 1 + + ix = ix(1) + (ix(2:end)-1)' * cumNr; + + end + + + % If the current closest point on the new grid has not been assigned yet, it is + % assigned the value of the current point on the old grid. Otherwise, the current + % closest point on the new grid is reassigned, if additionally, the distance to + % the current point on the old grid is smaller than the previously assigned one. + + if ~any(out_ix(1:i-1) == ix) + + out_s(ix) = s(i); + + if ~isempty(out_Ub) + + out_Ub(ix) = Ub(i); + + end + + out_ix(i) = ix; + + r_all(i) = r; + + else + + i_old = i_all(out_ix(1:i-1) == ix); + + if r < r_all(i_old) + + out_s(ix) = s(i); + + if ~isempty(out_Ub) + + out_Ub(ix) = Ub(i); + + end + + out_ix(i) = ix; + + r_all(i) = r; + + out_ix(i_old) = -out_ix(i_old); % eliminates from further comparisons + + else + + out_ix(i) = -out_ix(i); % eliminates from further comparisons + + end + + end + +end + + +end + diff --git a/MATLAB/library_m/f_ap_projected.m b/MATLAB/library_m/f_ap_projected.m new file mode 100644 index 0000000..181fcd8 --- /dev/null +++ b/MATLAB/library_m/f_ap_projected.m @@ -0,0 +1,326 @@ + +% C O P Y R I G H T N O T I C E +% +% Copyright ©2021. Institute of Science and Technology Austria (IST Austria). +% All Rights Reserved. The underlying technology is protected by PCT Patent +% Application No. PCT/EP2021/054650. +% +% This file is part of the AP demodulation library, which is free software: you can +% redistribute it and/or modify it under the terms of the GNU General Public License +% as published by the Free Software Foundation in version 2. +% +% This program 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 v2 for more details. You +% should have received a copy of the GNU General Public License v2 along with this +% program. If not, see https://www.gnu.org/licenses/. +% +% Contact the Technology Transfer Office, IST Austria, Am Campus 1, +% A-3400 Klosterneuburg, Austria, +43-(0)2243 9000, twist@ist.ac.at, for commercial +% licensing opportunities. +% +% All other inquiries should be directed to the author, Mantas Gabrielaitis, +% mgabriel@ist.ac.at + + + +function [m_out, e_out, iter] = f_ap_projected (s, Par, Ub, ix) +% > +% S H O R T D E S C R I P T I O N +% +% Calculates the modulator of a signal by using the AP-Basic algorithm. Signals +% defined in up to 3 dimensions are allowed. +% +% +% I N P U T A R G U M E N T S +% +% [s] - input signal. This is a real array of a chosen dimension D ≤ 3. +% +% [Par] - parameters characterizing the signal and demodulation procedure. This is a +% variable of the structure type. Its fields are as follows: +% +% .Fs - sampling frequencies for each dimension of the signal. This is an +% array of D elements. +% +% .Fc - cutoff frequencies of the modulator for each dimension of the signal. +% This is an array of D elements. +% +% .Et - infeasibility error tolerance used to control the termination of +% the demodulation algorithm. The iterative process is stopped when the +% infeasibility error, ϵ, drops to the level of .Et or below. If +% .Et ≤ 0, then the maximum allowed number of iterations, .Ni (see +% below), is completed. +% +% .Ni - maximum number of allowed iterations of the chosen AP algorithm. The +% chosen AP algorithm is iterated not more than .Ni times independent +% of whether ϵ drops at or below .Et. +% +% .Ns numbers of sample points of the original input signal (before any +% possible interpolation) in every dimension. This is an array with the +% number of elements equal to the number of signal dimensions. +% +% .Nx - numbers of sample points of the provided input signal (possibly +% interpolated) in every dimension. This is an array with the number +% of elements equal to the number of signal dimensions. +% +% .im - array with the iteration numbers at which the modulator estimates +% have to be saved for the output. If .im is empty, only the final +% modulator estimate is saved. The same is true when numel(i.im)=1 and +% .im(1)=.Ni. This field is optional (.im=[] is assumed by default). +% +% .ie - array with the iteration numbers at which the infeasibility error, ϵ, +% values have to be saved for the output. If .ie is empty, only the +% final error is saved. The same is true when numel(i.im)=1 and +% .im(1)=.Ni. This field is optional (.ie=[] is assumed by default). +% +% [Ub] - array with the upper bound constraints on the modulator. If Ub is an empty +% array, no upper bound constraint is assumed. Otherwise, Ub must be an array +% with the same dimensions as s. This input argument is optional (the default +% is Ub=[]). +% +% [ix] - linear indexes of the elements of an interpolated signal corresponding to +% every element of the original signal, i.e., mapping between the linear +% indexes of the original and interpolated signals. +% +% +% O U T P U T A R G U M E N T S +% +% [m_out] - modulator estimates at predefined interations (D+1 array). +% +% [e_out] - infeasibility error estimates at predefined iterations (1D array). +% +% [Ni_out] - number of iterations actually used (scalar). + + + +% [1] DEFINITION AND INITIALIZATION OF RELEVANT VARIABLES + +% Dimensionality of the input signal + +D = sum(size(s)~=1); + + + +% Numbers of sample points of the original and possibly interpolated signals + +ns = prod(Par.Ns); + +nx = prod(Par.Nx); + + + +% Full-wave rectification and scaling of the signal + +if D == 1 + + s = s(:); + +end + +s_abs = abs(s); % absolute value of the signal + +max_s = max(s_abs(:)); % maximum of the absolute-value signal + +s_abs = s_abs / max_s; % scaled absolute signal + +if ~isempty(Ub) + + Ub = Ub / max_s; + +end + + + +% Frequencies to filter out with the projection onto Mw + +fL = 1 + ceil(Par.Fc(:) ./ (Par.Fs(:)./Par.Nx)); + +fR = Par.Nx - fL + 2; + +i_cutoff = cell(D,1); + +for i=1:D + + i_cutoff{i} = fL(i)+1:fR(i)-1; + +end + + +% Initialization of the infeasibility error variable + +if Par.Et > 0 + + Etol = (Par.Et / max_s)^2; + +else + + Etol = Par.Et; + +end + +E = sum(s_abs(:).^2) / nx; + + + +% Initialization of the output arguments and related variables + +iter_m = 1; % iteration variable for storing the modulator estimates + +nim = numel(Par.im); + + +m_out = zeros([Par.Ns; max(numel(Par.im),1)]'); % output array with modulator estimates + +if ~isempty(Par.im) && Par.im(iter_m) == 1 + + m_out(1+(iter_m-1)*nx:iter_m*nx) = s_abs; + + iter_m = iter_m + 1; + +end + + + +iter_e = 1; % iteration variable for storing the infeasibility errors + +nie = numel(Par.ie); + + +e_out = zeros(numel(Par.ie),1); % output array with the infeasibility errors + +if ~isempty(Par.ie) && Par.ie(iter_e) == 1 + + e_out(iter_e) = E; + + iter_e = iter_e + 1; + +end + + + +% [3] CALCULATIONS + +iter = 0; + +m = s_abs; + +m_old = m; + +c = m; + + +while E > Etol && iter < Par.Ni + + iter = iter + 1; + + + % Projection onto Mw + + a = m; + + a = fftn(a); + + if D == 1 + + a(i_cutoff{1}) = 0; + + elseif D == 2 + + a(i_cutoff{1},:) = 0; + + a(:,i_cutoff{2}) = 0; + + elseif D == 3 + + a(i_cutoff{1},:,:) = 0; + + a(:,i_cutoff{2},:) = 0; + + a(:,:,i_cutoff{3}) = 0; + + end + + a = ifftn(a); + + + % Projection onto Cd + + m = a - c; + + m(m < s_abs) = s_abs(m < s_abs); + + if ~isempty(Ub) + + m(m > Ub) = Ub(m > Ub); + + end + + + % Infeasibility error + + E = (sum((m_old(:)-a(:)).^2) + sum((m(:)-a(:)).^2)) / (2*nx); + + + % Updated c and m_old + + c = c + m - a; + + m_old = m; + + + % Output (modulator) + + if iter_m <= nim && (iter == Par.im(iter_m) || (E <= Etol && nim == 1 && ... + Par.im(1) == Par.Ni)) + + m_out(1+(iter_m-1)*ns:iter_m*ns) = m(abs(ix)); + + iter_m = iter_m + 1; + + end + + + % Output (infeasibility error) + + if iter_e <= nie && (iter == Par.ie(iter_e) || (E <= Etol && nie == 1 && ... + Par.ie(1) == Par.Ni)) + + e_out(iter_e) = E; + + iter_e = iter_e + 1; + + end + +end + + +% [4] FINAL OUTPUT + +% If only the final estimate of the modulator is requested + +if isempty(Par.im) + + m_out(:) = m(abs(ix)); + +end + + +% Scaling back the modulator estimates + +m_out = max_s * m_out; + + +if isempty(Par.ie) + + + e_out = E; + +end + + +% Scaling back the error estimates + +e_out = max_s * sqrt(e_out); + + + diff --git a/MATLAB/library_mex/compilation.m b/MATLAB/library_mex/compilation.m new file mode 100644 index 0000000..edd6531 --- /dev/null +++ b/MATLAB/library_mex/compilation.m @@ -0,0 +1,66 @@ + + +% [1] Choose between 'dynamic' and 'static' linking. Currently, only the former is +% available + +CompType = 'dynamic'; + + + +% [2] Compilation flags + +CFlags = '''$CFLAGS -DAP_MEX_COMPILATION'''; + + + +% [3] Input and output filenames + +InFileName = 'f_ap_demodulation_mex.c'; + +OutFileName = './f_ap_demodulation_mex'; + + + +% [4] Path to the main directory of the Intel MKL library +% +% This path must be changed depending on where on your sustem the Intel MKL +% library is installed. + +IntelPath = '/usr/local/intel2019/compilers_and_libraries_2019/linux/mkl/'; + + + +% [5] Paths to the headers and executables of the AP demodulation and Intel MKL +% libraries + +IncludeDirAP = '/home/mgabriel/Documents/simulations/auditory_demodulation/code_for_sharing/C/library/'; + +IncludeDirMKL = [IntelPath, 'include']; + +LibraryDirMKL = [IntelPath, 'lib/intel64/']; + + + +% [6] external libraries to include for compilation + +Libraries = {[LibraryDirMKL 'libmkl_rt.so '], '-lm ','-ldl'}; + + + +% [7] Compilation + +if strcmpi(CompType,'Dynamic') + + Final = [ 'mex CFLAGS=', CFlags, ... + ' ', InFileName, ... + ' -I', IncludeDirAP, ... + ' -I', IncludeDirMKL, ... + ' ', Libraries{:}, ... + ' -output ', OutFileName ]; + +elseif strcmpi(CompType,'static') + % TODO +end + + +eval(Final) diff --git a/MATLAB/library_mex/f_ap_demodulation_mex.c b/MATLAB/library_mex/f_ap_demodulation_mex.c new file mode 100644 index 0000000..1f934f8 --- /dev/null +++ b/MATLAB/library_mex/f_ap_demodulation_mex.c @@ -0,0 +1,763 @@ + +/* C O P Y R I G H T N O T I C E + * + * Copyright ©2021. Institute of Science and Technology Austria (IST Austria). + * All Rights Reserved. The underlying technology is protected by PCT Patent + * Application No. PCT/EP2021/054650. + * + * This file is part of the AP demodulation library, which is free software: you can + * redistribute it and/or modify it under the terms of the GNU General Public License + * as published by the Free Software Foundation in version 2. + * + * This program 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 v2 for more details. You + * should have received a copy of the GNU General Public License v2 along with this + * program. If not, see https://www.gnu.org/licenses/. + * + * Contact the Technology Transfer Office, IST Austria, Am Campus 1, + * A-3400 Klosterneuburg, Austria, +43-(0)2243 9000, twist@ist.ac.at, for commercial + * licensing opportunities. + * + * All other inquiries should be directed to the author, Mantas Gabrielaitis, + * mgabriel@ist.ac.at + */ + + +/* C O N T E N T S + * + * A MEX gateway routine for f_ap_demodulation: + * + * (1) mexFunction. + */ + + + +#include + +#include "mex.h" + +#include "f_ap_demodulation.c" + + + +void mexFunction( int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) +{ +/* S H O R T D E S C R I P T I O N + * + * f_ap_demodulation_mex estimates the amplitude modulator of a given signal in the + * offline mode by using the AP demodulation approach formulated and developed in + * + * M. Gabrielaitis. "Fast and Accurate Amplitude Demodulation of Wideband Signals," + * IEEE Transactions on Signal Processing, vol. 69, pp. 4039-4054, 2021. DOI: + * 10.1109/TSP.2021.3087899. + * + * Please cite the above publication when referring to the algorithms implemented in + * this code or the code itself. + * + * See https://github.com/mgabriel-lt/ap-demodulation for the latest version of the + * code and user-friendly explanations on the working principle, domains of + * application, and advice on the usage of different AP demodulation algorithms in + * practice. + */ + +/* I N P U T A R G U M E N T S + * + * [1] - input signal. This is a real array with at least two elements and of a + * chosen dimension D ≤ 3. If the fourth input argument, t, is nonempty, i.e., + * interpolation of the input signal is required, then s must be a 1D array. + * + * [2] - parameters characterizing the signal and demodulation procedure. This is a + * variable of the structure type. Its fields are as follows: + * + * .Al - demodulation algorithm. Possible options are: 'B' - AP-Basic, + * 'A' - AP-Accelerated, 'P' - AP-Projected. + * + * .Fs - sampling frequencies for each dimension of the signal. This is an + * array of D elements. + * + * .Fc - cutoff frequencies of the modulator for each dimension of the signal. + * This is an array of D elements. + * + * .Et - infeasibility error tolerance used to control the termination of the + * demodulation algorithm. The iterative process is stopped when the + * infeasibility error, ϵ, drops to the level of .Et or below. If + * .Et ≤ 0, then the maximum allowed number of iterations, .Ni (see + * below), is completed. + * + * .Ni - maximum number of allowed iterations of the chosen AP algorithm. The + * chosen AP algorithm is iterated not more than .Ni times independent of + * whether ϵ drops at or below .Et. + * + * .Nr - numbers of sample points on the refined uniform interpolation grid in + * every dimension. This is an array of D elements. This field must be + * defined only if the fourth input argument, t, is provided (see below). + * Otherwise, .Nr is ignored. + * + * .Cp - compression parameter. If .Cp > 1, demodulation is performed by using + * signal compression (see equations (12)-(13) in the paper cited in the + * header of this file. This field is optional (the default is .Cp=1). + * + * .Br - indicator of premature termination of the 'AP-Accelerated' algorithm + * when the λ factor drops below one. If .Br=1 (recommended), premature + * termination is assumed (default). Otherwise, if .Br=0, the AP-A is not + * stopped prematurely even when λ decreases below 1. This field is + * required only if .Al='A'. It is optional (the default is .Br=1). + * + * .im - array with the iteration numbers at which the modulator estimates + * have to be saved for the output. If .im is empty, only the final + * modulator estimate is saved. This field is optional (.im=[] is assumed + * by default). + * + * .ie - array with the iteration numbers at which the infeasibility error, ϵ, + * values have to be saved for the output. If .ie is empty, only the + * final error is saved. This field is optional (.ie=[] is assumed by + * default). + * + * [3] - array with the upper bound constraints on the modulator. Two options are + * possible. If Ub is an empty array, no upper bound constraint is assumed. + * Otherwise, Ub must be an array with the same dimensions as s. This input + * argument is optional (Ub=[] is assumed by default). + * + * [4] - sampling coordinates of the input signal. This is a real 2D array with the + * number of columns equal to the dimension of the input signal and the number + * of rows equal to the number of signal sample points. This input argument is + * optional and should be used only if the input signal is sampled + * nonuniformly. If this input argument is provided, then the second input + * argument must contain the field .Nr (see above). + */ + +/* O U T P U T A R G U M E N T S + * + * [1] - modulator. + * + * [2] - infeasibility error of the modulator, if .Online == 0, or CPU time for each + * demodulation window, if .Online == 1. + * + * [3] - number of AP iterations used (the average number is returned, in the + * case of window splitting or online demodulation). + * + * [4] - running time of the program in seconds. + */ + + +/***********************************************************************************/ +/*********************** DECLARATIONS AND INITIALIZATIONS *************************/ +/***********************************************************************************/ + + + int exitflag = 0; + + mxChar *pr_c; + + double *pr_in; + + mxArray *pr_in2; + + mwSize D; + + const mwSize *Ns; + + const mwSize *sz_t; + + mwSize *Nm; + + long N; + + + long i; + + long n = 1; + + long n_im; + + long n_ie; + + double *pr_s; + + double *pr_t = NULL; + + double *Ub = NULL; + + double *e = NULL; + + double *pr_m; + + double *pr_e; + + long iter; + + + struct s_ParamsAP Par; + + + clock_t start, end; + + double cpu_time_used; + + + + + /* Start of measuring the running time */ + + start = clock(); + + + +/***********************************************************************************/ +/************************** READOUT OF INPUT ARGUMENTS ****************************/ +/***********************************************************************************/ + + + /* Checking if the number of input arguments is right */ + + if( nrhs < 2 || nrhs > 4 ) + + mexErrMsgIdAndTxt("AP_Demodulation:InpVal", "Two to four input arguments "\ + "are required!"); + + + /* Checking if the number of output arguments is right */ + + if( nlhs < 1 || nlhs > 4 ) + + mexErrMsgIdAndTxt("AP_Demodulation:InpVal", "One to four output arguments "\ + "are required!"); + + + + + /* Making sure that the 1st input argument is a real array with right dimensions */ + + N = mxGetNumberOfElements(prhs[0]); + + Ns = mxGetDimensions (prhs[0]); + + if ( nrhs < 4 || mxIsEmpty(prhs[3]) ) + { + D = mxGetNumberOfDimensions (prhs[0]); + + if ( D == 2 && (Ns[0] == 1 || Ns[1] == 1) ) + + D = D - 1; + } + else + { + sz_t = mxGetDimensions (prhs[3]); + + D = sz_t[1]; + } + + + if ( !mxIsNumeric(prhs[0]) || N < 2 || D > 3 ) /*what if complex ??? */ + + mexErrMsgIdAndTxt("AP_Demodulation:InpVal", "The 1st input argument must be"\ + " a signal with at least two sample points and no more"\ + " than three dimensions!"); + + + + + /* Getting a pointer to the 1st input argument */ + + pr_s = mxGetPr(prhs[0]); + + + + + /* Checking if the 2nd input argument is a structure */ + + if( !mxIsStruct(prhs[1]) ) + + mexErrMsgIdAndTxt("AP_Demodulation:InpVal", "The 2nd input argument must "\ + "be a structure!"); + + + + + /* Creating a local input parameter structure Par */ + + pr_in2 = mxGetField(prhs[1], 0, "Fs"); + + if ( pr_in2 == NULL ) + + mexErrMsgIdAndTxt("AP_Demodulation:InpVal", "The 2nd input argument must "\ + "have the field Fs!"); + + else if ( mxGetNumberOfElements(pr_in2) != D || !mxIsNumeric(pr_in2) ) + + mexErrMsgIdAndTxt("AP_Demodulation:InpVal", "Field 'Fs' of the 2nd input "\ + "argument must be a numeric array with the number of "\ + "elements equal to the dimension of the signal!"); + + else + + Par.Fs = mxGetPr(pr_in2); + + + + + pr_in2 = mxGetField(prhs[1], 0, "Fc"); + + if ( pr_in2 == NULL ) + + mexErrMsgIdAndTxt("AP_Demodulation:InpVal", "The 2nd input argument must "\ + "have the field Fc!"); + + else if ( mxGetNumberOfElements(pr_in2) != D || !mxIsNumeric(pr_in2) ) + + mexErrMsgIdAndTxt("AP_Demodulation:InpVal", "Field 'Fc' of the 2nd input "\ + "argument must be a numeric array with the number of "\ + "elements equal to the dimension of the signal!"); + + else + + Par.Fc = mxGetPr(pr_in2); + + + + + pr_in2 = mxGetField(prhs[1], 0, "Al"); + + if ( pr_in2 == NULL ) + + mexErrMsgIdAndTxt("AP_Demodulation:InpVal","The 2nd input argument must "\ + "have the field 'Al'!"); + + else if ( !mxIsScalar(pr_in2) || !mxIsChar(pr_in2) ) + + mexErrMsgIdAndTxt("AP_Demodulation:InpVal", "Field 'Al' of the 2nd input "\ + "argument must be a char scalar!"); + + else + { + pr_c = mxGetChars(pr_in2); + + Par.Al = pr_c[0]; + } + + + + + pr_in2 = mxGetField(prhs[1], 0, "Et"); + + if ( pr_in2 == NULL ) + + mexErrMsgIdAndTxt("AP_Demodulation:InpVal", "The 2nd input argument must "\ + "have the field Et!"); + + else if ( !mxIsScalar(pr_in2) || !mxIsNumeric(pr_in2) ) + + mexErrMsgIdAndTxt("AP_Demodulation:InpVal", "Field 'Et' of the 2nd input "\ + "argument must be a numeric scalar!"); + + else + + Par.Et = mxGetScalar(pr_in2); + + + + + pr_in2 = mxGetField(prhs[1], 0, "Ni"); + + if ( pr_in2 == NULL ) + + mexErrMsgIdAndTxt("AP_Demodulation:InpVal","The 2nd input argument must "\ + "have the field Ni!"); + + else if ( mxGetNumberOfElements(pr_in2) != 1 || !mxIsNumeric(pr_in2) ) + + mexErrMsgIdAndTxt("AP_Demodulation:InpVal", "Field 'Ni' of the 2nd input "\ + "argument must be a numeric scalar!"); + + else + + Par.Ni = (long) mxGetScalar(pr_in2); + + + + + pr_in2 = mxGetField(prhs[1], 0, "Nr"); + + if ( nrhs == 4 && !mxIsEmpty(prhs[3]) && pr_in2 == NULL ) + + mexErrMsgIdAndTxt("AP_Demodulation:InpVal", "The 2nd input argument must "\ + "have the field Nr!"); + + else if ( nrhs == 4 && !mxIsEmpty(prhs[3]) && \ + (mxGetNumberOfElements(pr_in2) != D || !mxIsNumeric(pr_in2)) ) + + mexErrMsgIdAndTxt("AP_Demodulation:InpVal", "Field 'Nr' of the 2nd input "\ + "argument must be a numeric array with the number of "\ + "elements equal to the dimension of the signal!"); + + else if (nrhs == 4) + { + pr_in = mxGetPr(pr_in2); + + Par.Nr = (long*) mxMalloc(D*sizeof(long)); + + if (Par.Nr == NULL) mexErrMsgIdAndTxt("AP_Demodulation:InpVal", \ + "Out of memory!"); + + for (i=0; i 2 && !mxIsEmpty(prhs[2]) ) + { + + n = (long) mxGetNumberOfElements(prhs[2]); + + + Ub = (double*) mxMalloc(N*sizeof(double)); + + if (Ub == NULL) mexErrMsgIdAndTxt("AP_Demodulation:InpVal",\ + "Out of memory!"); + + + if ( !mxIsNumeric(prhs[2]) || n != N ) + + mexErrMsgIdAndTxt("AP_Demodulation:InpVal", "The third input argument "\ + "must be either empty or a numeric array with the "\ + "same number of elements as the input signal!"); + + else + { + + pr_in = mxGetPr(prhs[2]); + + for (i=0; i 3 && !mxIsEmpty(prhs[3]) ) + { + + sz_t = mxGetDimensions (prhs[3]); + + if ( !mxIsNumeric(prhs[3]) || mxGetNumberOfDimensions(prhs[3]) != 2 || \ + sz_t[0] != N || sz_t[1] != D ) + + mexErrMsgIdAndTxt("AP_Demodulation:InpVal", "The 4th input argument "\ + "must be either empty or a 2D array with the number "\ + "of columns equal to the dimension of the input "\ + "signal and the number of rows equal to the number "\ + "of signal sample points!"); + + else + + pr_t = mxGetPr(prhs[3]); + + } + + + +/***********************************************************************************/ +/******************** MEMORY ALLOCATION FOR OUTPUT ARGUMENTS **********************/ +/***********************************************************************************/ + + + /* Allocating memory for the first two output arguments */ + + if ( nrhs < 4 || mxIsEmpty(prhs[3]) ) + { + Nm = (mwSize*) mxMalloc((D+1)*sizeof(mwSize)); + + if (Nm == NULL) mexErrMsgIdAndTxt("AP_Demodulation:InpVal","Out of memory!"); + + + for (i=0; i 1) + + plhs[1] = mxCreateDoubleMatrix(n_ie, 1, mxREAL); + + + + /* Getting pointers to the first two output arguments */ + + pr_m = mxGetPr(plhs[0]); + + if (nlhs > 1) + + pr_e = mxGetPr(plhs[1]); + + else + { + e = (double*) mxMalloc(n_ie*sizeof(double)); + + if (e == NULL) mexErrMsgIdAndTxt("AP_Demodulation:InpVal","Out of memory!"); + + + pr_e = e; + } + + + +/***********************************************************************************/ +/********************************* CALCULATION *************************************/ +/***********************************************************************************/ + + + /* Calling the computational routine */ + + exitflag = f_ap_demodulation (pr_s, &Par, Ub, pr_t, pr_m, pr_e, &iter); + + + +/***********************************************************************************/ +/************************** DYNAMIC MEMORY DEALLOCATION ****************************/ +/***********************************************************************************/ + + + mxFree(Par.Nr); + + mxFree(Par.im); + + mxFree(Par.ie); + + mxFree(Par.Ns); + + mxFree(Ub); + + mxFree(Nm); + + mxFree(e); + + + +/***********************************************************************************/ +/************************ READOUT OF THE COMPUTATION TIME **************************/ +/***********************************************************************************/ + + + /* Assigning the third output argument if requested */ + + if (nlhs > 2) + + plhs[2] = mxCreateDoubleScalar(iter); + + + + + /* Stopping measuring the running time */ + + end = clock(); + + cpu_time_used = ((double) (end - start)) / CLOCKS_PER_SEC; + + + + + /* Assigning the fourth output argument if requested */ + + if (nlhs > 3) + + plhs[3] = mxCreateDoubleScalar(cpu_time_used); + + + + if (exitflag != 0) mexErrMsgIdAndTxt("AP_Demodulation:InpVal", ""); + +} + diff --git a/README.md b/README.md new file mode 100644 index 0000000..f370cf0 --- /dev/null +++ b/README.md @@ -0,0 +1,13 @@ +# AP Demodulation +**AP Demodulation** is a mini-library of routines implementing the new approach to amplitude demodulation of arbitrary-bandwidth signals introduced in + +> [M. Gabrielaitis. "Fast and Accurate Amplitude Demodulation of Wideband signals," IEEE Transactions on Signal Processing, vol. 69, pp. 4039-4054, 2021. DOI: 10.1109/TSP.2021.3087899.](https://ieeexplore.ieee.org/document/9449954)[1](#footnote1) + +In this approach, demodulation is posed as a recovery problem of a signal from an unlabelled mix of its true and corrupted sample points. The recovery is achieved via tailor-made alternating projection algorithms. In the context of narrowband signals, the new method outperforms classical demodulation algorithms in terms of robustness to data distortions and compatibility with nonuniform sampling. When considering demodulation of wideband signals, it surpasses the current state-of-the-art techniques in terms of computational efficiency by up to many orders of magnitude. Such performance enables practical applications of amplitude demodulation in previously inaccessible settings. Specifically, online and large-scale offline demodulation of wideband signals, signals in higher dimensions, and poorly-sampled signals become practically feasible. + +#### Organization and usage of the library +Currently, the library includes C and MATLAB implementations, organized into separate folders in the repository root. The MATLAB implementation is divided into m-file and MEX function sections.[2](#footnote2) Folders of each programming language version contain comprehensive README.txt files that describe the contents of these folders and the usage and compilation guidelines of the library. In all cases, code with different signal demodulation examples is included to illustrate application of the library in practice. Note that the functionality of each library interface is accessible by a user via only a single function. +___ +1: A freely available preprint of this work is on [arXiv](https://arxiv.org/abs/2102.04832). + +2: The MEX implementation features shorter computing times than its m-file counterpart. The difference between the two versions is most noticeable for short 1D signals, when it may reach ~100 times. For very long signals, the MEX version advantage saturates at ~4 times. For signals defined in higher dimensions, the difference between the two implementations is small.