Fast and Accurate Amplitude Demodulation of Wideband Signals
11 Amplitude Demodulation of Wideband Signals
Mantas Gabrielaitis
Abstract —Amplitude demodulation is a classical operationused in signal processing. For a long time, its effective applicationsin practice have been limited to narrowband signals. In this work,we generalize amplitude demodulation to wideband signals. Wepose demodulation as a recovery problem of an oversampled cor-rupted signal and introduce special iterative schemes belonging tothe family of alternating projection algorithms to solve it. Sensiblychosen structural assumptions on the demodulation outputs allowus to reveal the high inferential accuracy of the method over arich set of relevant signals. This new approach surpasses currentstate-of-the-art demodulation techniques apt to wideband signalsin computational efficiency by up to many orders of magnitudewith no sacrifice in quality. Such performance opens the doorfor applications of the amplitude demodulation procedure in newcontexts. In particular, the new method makes online and large-scale offline data processing feasible, including the calculation ofmodulator-carrier pairs in higher dimensions and poor samplingconditions, independent of the signal bandwidth. We illustratethe utility and specifics of applications of the new method inpractice by using synthetic and natural speech signals.
Index Terms —Alternating projections, amplitude demodula-tion, convex programming, fast algorithms, multidimensionalsignals, speech processing, wideband signals.
I. I
NTRODUCTION A MPLITUDE demodulation refers to the decompositionof a signal into a product of a slow-varying modulator-envelope and a fast-varying carrier. First introduced in radiocommunications [1], this procedure has found applications indata acquisition and processing related to a broad range ofphenomena. Automatic speech recognition [2], atomic forcemicroscopy [3], ultrasound imaging [4], brainwave [5], seismictrace [6], and fingerprint [7] analyses are a few among manyexamples to mention.Originally, amplitude demodulation was intended for usewith signals built of locally sinusoidal, i.e., narrowband, carri-ers. Several classical approaches excel in this setting, with Ga-bor’s analytic-signal (AS) method being a long-standing cham-pion [8], [9]. Nonetheless, many relevant problems inevitablyrequire demodulating signals that feature wideband carriers,typically of (quasi)-harmonic, (quasi)-random, or spike-trainorigin [10]–[18] (see Suppl. Mat. H for an overview). Whenapplied to them, the classical techniques fail, misleadinglymixing the carrier and modulator information [19], [20].For a long time, no consistent and accurate way of de-modulating wideband signals was known. Typically, a proxyof the modulator would be obtained by rectifying and thenlow-pass filtering the signal. Different implementations ofthis procedure, each adapted for a specific signal class, weresuggested (see, e.g., [17], [21], [22]). The estimates of signalmodulators obtained in this way, however, are neither accurate
M. Gabrielaitis is with the Institute of Science and Technology Austria,3400 Klosterneuburg, Austria (e-mail: [email protected]) nor consistent between different methods. The carriers andmodulators are not appropriately separated either, i.e., they canbe demodulated further by iterating the same procedure [23].Moreover, the carrier estimates are often unbounded, even inwell-defined situations (see, e.g., [23, Fig. 3.1]).Recently, two powerful demodulation approaches suitableto signals with arbitrary bandwidths have been formulated.Turner and Sahani shaped demodulation into a statisticalinference problem [23], [24]. In this so-called probabilisticamplitude demodulation (PAD) approach, the modulator andcarrier are inferred from the signal as latent variables of anappropriately selected statistical model. Mathematically, PADdefines a maximization of a posteriori probability, a high-dimensional nonlinear optimization task. In another work, Selland Slaney chose a deterministic route to demodulation [19].In their linear-domain convex (LDC) approach, the modulatoris described as a minimum-power signal with penalized high-frequency terms lying above the original waveform. This prob-lem is convex and thus amenable to more efficient optimizationmethods than the PAD.The PAD and LDC techniques separate the modulator andcarrier information of various synthetic wideband signals witha high degree of accuracy [19], [23]. The principal weaknessof these approaches is a huge associated computational burden,which impedes their use in practical situations (see Section IVfor the performance evaluations). In particular, online or large-scale offline signal processing is out of reach for the PAD andLDC demodulations. Besides, derivations of these methods areguided more by high-level modulator or carrier properties andcomputational tractability rather than strict recovery condi-tions. Hence, the boundaries of their validity in the contextof real-world signals are somewhat blurred.In this work, we frame demodulation as a problem of modu-lator recovery from an unlabeled mix of its true and corruptedsample points. We show that, under some loose constraints oncarriers and modulators, high-accuracy demodulation can beachieved through exact or approximate norm minimization.We introduce different versions of custom-made alternatingprojection algorithms and test them in numerical experimentsto solve this task. The new approach is shown to be free ofthe performance limitations inherent to the PAD and LDCmethods. In particular, it combines the computational economyof the classical AS technique with a capacity to recover a widerange of arbitrary-bandwidth signals. We reveal the power ofthe new approach in terms of efficiency, accuracy, consistency,and robustness to corrupted data through theoretical analysisand illustrate it using synthetic signals with known structure.The use of the new method in realistic online and offlinesettings is demonstrated by applying it to natural speech. a r X i v : . [ ee ss . SP ] F e b II. M
ATHEMATICAL F ORMULATION OF THE P ROBLEM
In what follows, we assume the representation of a real-valued signal s ( t ) formed by a finite collection of its valuesuniformly sampled over a limited time interval: s i ≡ s ( t i ) , i ∈ I n = { , , . . . , n } . Thus, a realization of the signal, s ≡ ( s , s , . . . , s n ) T , is an element of an n -dimensional Eu-clidean space R n , i.e., a linear space equipped with an innerproduct h s (1) , s (2) i = P ni =1 ( s (1) i · s (2) i ) , which induces theEuclidean norm k s k = p h s , s i . A. Demodulation constraints
The task of demodulation is to factorize a signal s ∈ R n into a modulator m ∈ R n and a carrier c ∈ R n : s = m ◦ c , (1)where symbol ◦ denotes an elementwise product of twovectors. There exists an uncountable number of pairs of m and c that satisfy (1). Thus, further constraints are needed todefine its unique solution. It is precisely these constraints thatgive a distinct character to different demodulation methodsand set the domain of their validity [9], [19], [23], [25].In this work, we introduce the extra demodulation restric-tions by imposing some general assumptions on m and c .We define feasible modulators as elements of a convex set M ω = S ≥ ∩ S ω , (2)where S ≥ = { x ∈ R n : x i ≥ , i ∈ I n } , S ω = { x ∈ R n : ( Fx ) i = 0 , i ∈ ( I n \ I ωn ) } , (3) I ωn = { i ∈ I n : i ≤ ω } ∪ { i ∈ I n : i > n + 1 − ω } . In (3), F denotes the operator of the unitary discrete Fouriertransform (DFT), and ( . . . ) i marks the i -th component of theargument vector. Hence, in our framework, modulators arenonnegative low-pass signals whose rate of variation is limitedby the cutoff frequency ≤ ω ≤ d n/ e . This is a formaldefinition of the classical modulator-envelope [1], [26].We declare feasible carriers as elements of a nonconvex set C d = S | .. |≤ ∩ S { } ,d , (4)where S | .. |≤ = { x ∈ R n : | x i | ≤ , i ∈ I n } , S { } ,d = (cid:8) x ∈ R n : ( ∀ i ) P i + d − j = i I { } ( | x j | ) ≥ , ( ∃ i ) P i + d − j = i I { } ( | x j | ) = 1 (cid:9) , (5)with I { } being the indicator function of the singleton { } .The set S | .. |≤ implies boundedness of c between − and .This restriction follows from the standard notion that the time-dependent amplitude of an amplitude-modulated s is purely setby m . Meanwhile, S { } ,d fixes to d the maximum gap betweenany two neighboring components of c whose absolute valuesare equal to . Note that C d ∩ C d = ∅ if d = d . As shownnext, this constraint allows formulating extensive demodu-lation guarantees while only moderately affecting the scopeof relevant carriers. Bandwidth-wise, C d covers the wholerange, from zero (sinusoidal) to flat (random spike) bandwidthsignals, and defines the qualifier “wideband” used in this work. B. Demodulation as modulator recovery
Note that, assuming c ∈ C d , | s | can be seen as a corruptedversion of m : | s i | = m i when | c i | = 1 , and | s i | 6 = m i otherwise. Further, if m can be found from | s | , c followsfrom (1) uniquely, except the sample points with m i = 0 . Thelatter, if any, are sparse and can be typically interpolated fromthe neighboring points. Hence, in our case, demodulation isvirtually a problem of reconstructing m from a mix of its true( | s i | = m i ) and corrupted ( | s i | 6 = m i ) sample points when theclass of each point is unknown. This viewpoint is at the coreof the developments that follow next. C. Modulator recovery through norm minimization
Our approach to demodulation builds around the estimator ˆm = arg min x ∈S ≥| s | ∩S $ k x k , (6)where S ≥| s | = { x ∈ R n : x i ≥ | s i | , i ∈ I n } . Note that m ∈S ≥| s | ∩ S $ if $ ≥ ω . The restriction S ≥| s | assures that x doesnot fall below the true sample points of m . If, besides, thetrue sample points are spread densely enough, we expect thenorm minimization to enforce ˆ m i = m i at these points. Butthen, ˆm = m by the discrete sampling theorem. A foundationfor this intuitive consideration is laid by the following results. Proposition II.1 (proof in Suppl. Mat. B) . For almost every m ∈ M ω , ˆm = m only if c ∈ C d , $ ≥ ω , and n s ≡ P ni =1 I { } ( | c i | ) ≥ $ + ω − . Proposition II.2 (proof in Suppl. Mat. B) . Consider m ∈ M ω and ˜c ∈ C ˜ d with | ˜ c i | = 1 for i ∈ J n ⊆ I n . If ˆm = m holdsfor the m and ˜c , then it also holds for every pair made of thesame m and any c ∈ C d with | c i | = 1 for i ∈ J n . Proposition II.3 (proof in Suppl. Mat. B) . Assume m ∈ M ω and c ∈ C d with $ ≥ ω . If, additionally, there exist d ∈ I n and i ∈ I d such that n s ≡ ( n/d ) ∈ N + , n s ≥ $ + ω − , and | c i +( j − · d | = 1 for every j ∈ I n s , then ˆm = m .Proposition II.1 reveals a tight match of ˆm to C d : no m ∈M ω can be inferred from s by ˆm precisely if c / ∈ C d . Italso establishes the central role of the presence of true samplepoints in the recovery: for almost every m ∈ M ω , at least $ + ω − such points are needed. Proposition II.2 furtherconsolidates the latter view by stating that the success of theexact recovery of an m ∈ M ω via ˆm depends only on thenumber and positions of the true sample points.In Proposition II.1 , $ ≥ ω and n s ≥ $ + ω − imply n s ≥ ω − , which is a sufficient condition for m recovery in the classical setup when all true sample points areknown (see the remark below the proof of Proposition A.1 in Suppl. Mat. A). Hence, the data corruption manifesting inour problem necessitates further constraints on the number orpositions of true sample points. In particular,
Proposition II.3 certifies a full recovery of m if $ ≥ ω , and there exists a (notnecessarily known) subset of at least $ + ω − regularly-spacedtrue sample points. The latter condition covers a wide rangeof practically relevant carriers, including: (1) the classical sin(2 πν t + φ ) with ν ≥ ω , (2) harmonic signals, (3) regular spike-trains of | c i | = 1 . More generally, any (non)stationarytime-series with regularly placed | c i | = 1 regardless of theremaining points are eligible.In addition to the regularity of true sample points, Proposi-tion II.3 requires n/d to be an integer. Nevertheless, numericalexperiments reveal that both of these conditions can be ignoredwithout practically relevant consequences (see Suppl. Mat. Cand Fig. 10 there). In particular, we found that the discrep-ancy between m and ˆm is vanishing with an overwhelmingprobability for any c ∈ C d if $ ≥ ω , and d n/d e ≥ $ − .This result noticeably extends the scope of recovery conditionsover the domain of practically relevant (quasi-)regular andstochastic carriers. Among the examples are nonstationarysinusoidal and harmonic signals and arbitrary spike-trains withthe distance between neighboring spikes at or below d points.Note that n s ≥ d n/d e . Hence the relaxation of the strictregularity condition on | c i | = 1 sample points comes atthe expense of a slightly tighter constraint on n s : compare n s ≥ $ − vs. n s ≥ $ + ω − . The latter statement isexact and is established as an intermediate result in the proofof Proposition II.1 .Another important generalization of the recovery conditionscomes with the following inequality:
Proposition II.4 (proof in Suppl. Mat. B) . Consider m ∈ M ω and c ∈ S | .. |≤ . Take n s ≥ $ − sample points of s = m ◦ c whose indexes are defined as entries of any chosen r ∈ N n s + with r i +1 − r i = n/n s for every i ∈ I n s . Then, k m − ˆm k / k m k ≤ q − P n s i =1 s r i / P n s i =1 m r i . (7)Hence, if one can find a sequence of at least $ − regularly-spaced sample points with | s i | sufficiently close to m i , thenthe relative recovery error is close to 0 in terms of (7). Thisresult endows ˆm with the stability to discrepancies from therecovery conditions discussed earlier. At the same time, itprovides approximate recovery guarantees for a wider range ofstochastic and (quasi-)regular carriers besides those with fairlydensely-packed | c i | = 1 sample points. Due to the low-passrestriction on m , (7) is expected to hold approximately for anirregular r ∈ N n s + with r i +1 − r i ≤ d n/n s e as well.We finally note that, whereas ω and d are fixed propertiesof m and c , $ is a control parameter that must be specified.An appropriate $ , which satisfies the recovery conditions for-mulated above, can only be selected by using prior knowledgeon m and c or found in a supervised learning setup. D. Relaxation of the exact minimum-norm requirement
The norm-minimizing property of ˆm in (6) is criticalin formulating sharp recovery conditions. However, from apractical point of view, little would be lost if another estimator ˆm with only slightly larger than the minimum norm amongall elements of S ≥| s | ∩ S $ is used. Thus, we relax (6) tofind ˆm ∈ S ≥| s | ∩ S $ subject to k ˆm k ’ arg min x ∈S ≥| s | ∩S $ k x k (8)To specify the otherwise ambiguous relation operator ’ , werequest that ˆm obtained through (8) recovers m exactly, i.e., is norm-minimizing, for sinusoidal, harmonic, and spike-traincarriers covered by Proposition II.3 . As we see later, thisrestriction regularizes the numerical algorithms formulated inthe present work for sufficiently accurate demodulation wellbeyond those three classes of c . E. Method of solution
The algorithms that we introduce to solve (6) and (8) in thiswork fall in the domain of the so-called methods of alternatingprojections (APs). The defining feature of each AP method isan iterative calculation of a feasible point ( ˆm ∈ S ≥| s | ∩ S $ inour case) through alternating metric projections of its currentestimate onto the separate constraint sets. Initially proposedby von Neumann for two closed subspaces [27], this approachwas later extended to arbitrary closed convex sets of a Hilbertspace (see [28] for a review). Various implementations of theAP algorithms exist, featuring different domains of applica-tion, rate of convergence, and additional requirements satisfiedby the solutions [28], [29].We provide a rigorous mathematical basis on which theAP algorithms for solving the demodulation problem relyin Suppl. Mat. D, E, F. For a practical comprehension of thematerial that follows next, it is sufficient to know that: • The sets S ≥| s | and S $ are closed and convex. • A metric projection, or simply a projection henceforth,of z ∈ R n onto a closed convex set S ⊂ R n is a unique x z ∈ S with the smallest distance, i.e., k x z − z k , from z . • The projections onto S ≥| s | and S $ are respectivelyachieved by operators P S ≥| s | [ z ] = | s | + ( z − | s | ) ◦ θ ( z − | s | ) (9)and P S $ [ z ] = ( F − W $ F ) z . (10)Here, θ ( . . . ) is the Heaviside step function evaluated ele-mentwise. W $ is a diagonal matrix such that ( W $ ) ii =1 if i ∈ I $n , and ( W $ ) ii = 0 otherwise.To emphasize the nature of the underlying numerical algo-rithms, we name our new approach as AP demodulation. F. Relation to other problems and approaches
Demodulation is a counterpart of a more widely knownand studied problem of blind deconvolution: s = m ~ c .Indeed, both tasks admit the algebraic form of each other inthe Fourier domain. Nevertheless, the properties of m and c inherent to practical instantiations of amplitude demodulationand blind deconvolution differ significantly. These differencespredetermine the need for distinctive strategies to solve therespective tasks, as discussed next.One of the most powerful convex-programming-based de-convolution approaches, introduced in [30], builds on theassumption that m and c belong to known low-dimensionalsubspaces. There, recovery of m and c is achieved by minimiz-ing the nuclear, atomic, ‘ , or ‘ , norms of their outer product(in the subspace representation) subject to linear measurementconstraints of s [30]–[33]. This scheme successfully solvesmany practically relevant blind deconvolution cases, such as t (s) x mm (1) a (1) m (0) t (s) x mm (3) a (3) m (2) t (s) x mm (2) a (2) m (1) Fig. 1. The first three iterations of the AP-B algorithm applied to an amplitude-modulated sinusoidal signal. m stands for the real modulator. image deblurring, multipath channel protection, and super-resolution microscopy [30], [33]. However, the low-dimensionsubspace assumption, a crucial prerequisite of the approach,renders it inapt to deal with realistic carriers in the amplitudedemodulation context. Indeed, even a sinusoidal carrier witha fluctuating phase is hardly representable in this frame, notto mention more complex wideband signals met in practice.Moreover, the subspace model of m and c does not allowenforcing the amplitude contents to m exclusively.Deconvolution problems have also been approached byusing AP-like methods [34]–[37]. A general strategy of theexisting algorithms is to achieve deconvolution by an iterativerefinement of both m and c upon the requirement of exact[34], [36] or approximate [35], [37] adherence to the definingequality s = m ~ c and the support region, intensity range, andspectrum constraints implied on m and c or r = s − m ~ c .These methods differ significantly between themselves. Eachof them achieves satisfactory recovery by a judicious com-bination of specific constraint sets and the iterative schemeadjusted to specific classes of m and c . The nonconvexity of C d and the absence of efficient explicit projections onto thisset makes the application of the known deconvolution methodsunsuitable to amplitude demodulation. None of the currentAP-like deconvolution methods allow assigning the amplitudecontents to m purely either.Our formulation of the amplitude demodulation problemin Section II-B reveals it as a generalization of the classicaltask of band-limited signal recovery from true sample points.An AP method known under the name Papoulis-Gerchbergand its variants have been successfully applied in the lattersetting (see [38] for a review). The differences in the availableinformation on the recoverable signal lead to distinct strategiesin algorithmic approaches to these two problems. In particular,the Papoulis-Gerchberg methods rely entirely on known truedata. Thus, they are impossible to use for demodulationpurposes. The AP algorithms introduced in the present workcan be applied in the classical setting. However, not using theavailable information about the true data makes them inferiorto their classical counterparts, unless the sample points arefairly uniformly spread, as discussed in Section II-C.The approach suggested in the present work also has someparallels with the LDC demodulation method by [19]. Inthe latter case, (1) is accompanied by a constraint on themodulator m expressed as the solution of a convex quadraticprogramming problemminimize k w ◦ Fm k + k m k subject to | s i | ≤ m i ≤ max[ s ] ∀ i ∈ I n , (11) where w denotes the weighting vector. (11) was introducedheuristically, trying to quantify the intuitive notion of themodulator-envelope as a signal wrapping s from above.Practical applications suggest the LDC method defined by(1) and (11) being computationally most efficient and preciseamong all current techniques designed for demodulating sig-nals unreachable to classical algorithms [19], [24]. Thus, weuse it as a reference when evaluating the performance of thenewly-formulated approach of the present work.III. D EMODULATION A LGORITHMS
In this section, we formulate three algorithms representingthe core arsenal of the AP approach to demodulation. Sim-plicity, efficiency, and estimation accuracy of the algorithmsare the main aspects under consideration. We refer the readerto Suppl. Mat. F for proofs of all propositions found here.
A. AP-Basic
We start with the simplest possible AP algorithm, thereforenamed “AP-Basic” (AP-B).
Algorithm:
AP-Basic (AP-B) Set: N iter , (cid:15) tol Initialize: i = 0 , (cid:15) (0) = k s k / √ n , m (0) = | s | , a (0) = while (cid:15) ( i ) > (cid:15) tol and i < N iter do i = i + 1 a ( i ) = P S $ [ m ( i − ] m ( i ) = P S ≥| s | [ a ( i ) ] (cid:15) ( i ) = k m ( i ) − a ( i ) k / √ n end while Finalize: ˆm = m ( i ) Here, N iter stands for the maximum number of algorithmiterations. (cid:15) ( i ) is the infeasibility error at the i -th iteration,which is used to control the termination of the algorithm.Specifically, (cid:15) ( i ) measures the distance of the estimate of themodulator m ( i ) ∈ S ≥| s | from S $ and sets a lower boundon the convergence error: (cid:15) ( i ) ≤ k m ( i − − m † k / √ n (seeSuppl. Mat. G). The iterative process is stopped when (cid:15) ( i ) drops to the level of a predetermined threshold (cid:15) tol > orbelow. (cid:15) tol ≤ would force the completion of all N iter iterations of the algorithm. ˆm denotes the final estimate of themodulator. ˆm arbitrarily close to S ≥| s | ∩ S $ can be reached if N iter is sufficiently large: Proposition III.1.
A sequence m (0) , m (1) , . . . , m ( i ) , . . . formed by the AP-B algorithm for (cid:15) tol = 0 and N iter → + ∞ converges to some m † ∈ S ≥| s | ∩ S $ . The convergence is geometric and monotonic, i.e., there exist γ > and < r < such that k m ( i ) − m † k ≤ γ · r i and k m ( i +1) − m † k ≤k m ( i ) − m † k for i ≥ . It can be shown by example that the AP-B does not alwaysprovide minimum-norm estimators ˆm . However, it is expectedto do so at least approximately if some conditions are met. Weclarify this next with the help of Fig. 1, which displays the firstthree iterations of the AP-B applied to an example signal.First, note that the starting point m (0) is elementwise not-higher than the real modulator m (black). P S $ maps m (0) to a (1) , which, by definition of a metric projection, is its bestmean-squared-error (MSE) approximation in S $ (blue). By thedefinition of S $ , a (1) is nearly constant over time windowsshorter than n/ (2 π$ ) points. In general, the best constantMSE estimator of a sample of numbers is its average. Thus,as the best MSE estimator of m (0) over S $ , a (1) approximatesthe local average of m (0) values in a window of ≈ n/ (2 π$ ) points at every moment. If c ∈ S | .. |≤ , $ ≥ ω , and ≈ n/ (2 π$ ) sample points are sufficient to average out local variations of c , a (1) is supposed to be proportional to m , at least roughly.The first iteration is completed by the projection of a (1) backonto S ≥| s | to obtain m (1) (red).Applying the same reasoning as above, we deduce that, witheach iteration, a ( i ) , and thus m ( i ) , approaches m elementwise(see Fig. 1). In general, m ( i ) may exceed the level of the realmodulator m over time windows longer than ≥ n/ (2 π$ ) points for higher i before m † ∈ S ≥| s | ∩S $ is reached. However,as follows from the considerations of the previous paragraph,such segments of m ( i ) would be approximately compatiblewith S ≥| s | ∩ S $ and would not be considerably affected insubsequent iterations. Hence, ˆm obtained by the AP-B isexpected to follow the true sample points of m tightly. If thenumber of these points is sufficient, then ˆm ≈ m as well.The basis for the above considerations is laid by the factthat they are exact for some important types of carriers: Proposition III.2.
Consider m ∈ M ω and c ∈ C d with | c j | = P n/νk =1 (˜ c ν · k · e ı πν ( k − j − /n ) , where ˜ c ν · k ∈ C and n/ν ∈ N . If $ ≥ ω and ν ≥ $ + ω − , then a sequence m (0) , m (1) , . . . , m ( i ) , . . . formed by the AP-B algorithm for (cid:15) tol = 0 and N iter → + ∞ converges to m . Among others,
Proposition III.2 encompasses the sinusoidal,harmonic, and regular spike-train carriers covered by
Propo-sition II.3 . Thus, in these cases, AP-B satisfies the minimum-norm property, i.e., provides ˆm converging to a solution of(8). The condition ν ≥ $ + ω − in Proposition III.2 playsthe role of the inequality n/d ≥ $ + ω − in Proposition II.3 . B. AP-Accelerated
One of the potential weak points of the AP algorithms basedon pure alternating projections onto convex sets, like the AP-B,is relatively slow convergence [39]–[41]. Indeed, despite thegeometric nature of the convergence, the actual number ofiterations necessary to reach a specific error level may bearbitrarily large if the factor r in k m ( i ) − m † k ≤ γ · r i issufficiently close to 1. To address this issue, various acceler-ation schemes of the AP algorithms have been suggested for specific classes of the constraint sets [39], [42], [43]. Here,we propose a parameter-free accelerated version of the AP-Balgorithm specifically adapted for the demodulation problem.We refer to it as “AP-Accelerated” (AP-A). Algorithm:
AP-Accelerated (AP-A) Set: N iter , (cid:15) tol Initialize: i = 0 , (cid:15) (0) = k s k / √ n , m (0) = | s | , a (0) = while (cid:15) ( i ) > (cid:15) tol and i < N iter do i = i + 1 b ( i ) = P S $ [ m ( i − − a ( i − ] λ = k m ( i − − a ( i − k / k b ( i ) k a ( i ) = a ( i − + λ · b ( i ) m ( i ) = P S ≥| s | [ a ( i ) ] (cid:15) ( i ) = k m ( i ) − a ( i ) k / √ n end while Finalize: ˆm = m ( i ) Proposition III.3.
A sequence m (0) , m (1) , . . . , m ( i ) , . . . formed by the AP-A algorithm for (cid:15) tol = 0 and N iter → + ∞ converges to some m † ∈ S ≥| s | ∩ S $ . The convergence ismonotonic, i.e., k m ( i +1) − m † k ≤ k m ( i ) − m † k for i ≥ . Note that λ > except when P S $ is the identity operator,i.e., the trivial case of m = | s | . Indeed, it follows from thedefinition of P S $ [see (10)] and the unitary property of F that k m ( i − − a ( i − k > k P S $ [ m ( i − − a ( i − ] k = k b ( i ) k , if P S $ is not the identity operator. It is easy to see that ( a ( i ) − a ( i − ) = λ · ( P S $ [ m ( i − ] − a ( i − ) in the above algorithm.Moreover, if λ is fixed to 1 by force, the AP-A and AP-Balgorithms become identical. Therefore, the AP-A producesincrements from a ( i − to a ( i ) that are scaled up comparedwith those that were obtained by applying the AP-B algorithmfor the same iterations.To understand the working principle of the AP-A better,recall that P S $ [ m ( i ) ] , and thus a ( i ) , are nearly constantover time windows consisting of < n/ (2 π$ ) points (seeSection III-A). For a semiquantitative analysis, we can as-sume that this holds exactly. Let us denote a segment of ( m ( i − − a ( i − ) restricted to such a window by z . Then, b ( i ) defined in the same window is just ( l − · P lj =1 z j ) · ,and k m ( i − − a ( i − k corresponds to P lj =1 z j , where, l = b n/ (2 π$ ) c . Consequently, λ · b ( i ) , i.e., the incrementfrom a ( i − to a ( i ) , is given by (cid:0) P lj =1 z j / P lj =1 z j (cid:1) · . Itfollows from m ( i − = P S ≥| s | [ a ( i − ] that z is elementwisenonnegative. Therefore, (cid:0) P lj =1 z j / P lj =1 z j (cid:1) ≤ max[ z ] .However, max[ z ] corresponds to the difference between thereal modulator and a ( i − in the considered time window atleast approximately, if d n/d e ≥ $ − . Thus, while up-scaling a ( i ) − a ( i − at each iteration to accelerate the convergence,the AP-A also ensures that a ( i ) stays approximately within thebounds of the real modulator m . This property ensures that ˆm tightly follows m if the same conditions as required by theAP-B are met.We further note that (cid:0) P lj =1 z j / P lj =1 z j (cid:1) = max[ z ] , i.e., a ( i ) reaches m in a single iteration, if all but one element of z are equal to zero. Importantly, approximately this situation is faced in practice, as illustrated in Fig. 1. Specifically, withincreased i , ( m ( i − − a ( i − ) becomes mainly flat withonly a few separate elements considerably above 0 over timewindows shorter than n/ (2 π$ ) points. For comparison, theanalogous increment from a ( i − to a ( i ) is moderate andequals max[ z ] /l in the case of the AP-B method. Theseconsiderations explain the substantial speed-up provided bythe AP-A algorithm in practice. They also reveal that anyadditional acceleration steps in the AP-A would result inoverscaled ˆm , hence reducing the demodulation quality.The AP-A algorithm repeats the AP-B in terms of exactrecovery guarantees of Proposition III.2 : Proposition III.4.
Consider m ∈ M ω and c ∈ C d with | c j | = P n/νk =1 (˜ c ν · k · e ı πν ( k − j − /n ) , where ˜ c ν · k ∈ C and n/ν ∈ N . If $ ≥ ω and ν ≥ $ + ω − , then a sequence m (0) , m (1) , . . . , m ( i ) , . . . formed by the AP-A algorithm for (cid:15) tol = 0 and N iter → + ∞ converges to m . This result substantiates the semiquantitative argumentation ofthe AP-A convergence properties provided above and estab-lishes the respective ˆm as a numerical solution of (8). C. AP-Projected
As argued above, the AP-A and AP-B algorithms producemodulator estimates that are expected to tightly follow theoriginal m if the conditions analogous to those discussedin Section II-C are met. These estimates, however, do notnecessarily hold the minimum-norm property (6). A classicalAP scheme that allows achieving the minimum-norm solutionsis known under the name of Dykstra [44], [45]. In particular,Dykstra’s algorithm calculates the projection of a point in R n onto the feasible set. Thus, by choosing an appropriate initialcondition, the solution with a minimized norm can be obtained(see Proposition III.5 next and its proof in Suppl. Mat. F).We consider a version of this algorithm adapted to solve thedemodulation problem and call it “AP-Projected” (AP-P).
Algorithm:
AP-Projected (AP-P) Set: N iter , (cid:15) tol Initialize: i = 0 , (cid:15) (0) = k s k / √ n , m (0) = c (0) = | s | while (cid:15) ( i ) > (cid:15) tol and i < N iter do i = i + 1 a ( i ) = P S $ [ m ( i − ] m ( i ) = P S ≥| s | [ a ( i ) − c ( i − ] c ( i ) = m ( i ) − ( a ( i ) − c ( i − ) (cid:15) ( i ) = p ( k m ( i − − a ( i ) k + k m ( i ) − a ( i ) k ) / (2 · n ) end while Finalize: ˆm = m ( i ) Proposition III.5.
A sequence m (0) , m (1) , . . . , m ( i ) , . . . formed by the AP-P algorithm for (cid:15) tol = 0 and N iter → + ∞ converges to a unique m † ∈ S ≥| s | ∩ S $ such that k m † k ≤k x k for every x ∈ S ≥| s | ∩ S $ . The convergence is monotonic,i.e., k m ( i +1) − m † k ≤ k m ( i ) − m † k for i ≥ . The AP-P differs from the AP-B in that, before projecting apoint onto S ≥| s | , an increment produced by the projection onto this set in the previous iteration is subtracted. This correctionmay cause the infeasibility error k m ( i ) − a ( i ) k / √ n estimatedafter projecting onto S ≥| s | to drop to zero intermittently beforethe final solution is reached, making it an inappropriate optionas the stopping criterion. Hence, in contrast to the AP-B andAP-A algorithms, we defined the (cid:15) for the AP-P as a combina-tion of the infeasibility errors evaluated after projecting ontoboth sets S $ and S ≥| s | at each iteration (see line 8 above). Thiserror measure is strictly positive and converges to zero when N iter → + ∞ [46].The understanding of the convergence rate of Dykstra’sscheme is limited. It was shown that the convergence isgeometric for an intersection of half-spaces [47], [48]. Nev-ertheless, no equivalent result exists for other convex sets.Moreover, it was demonstrated that the convergence rate ofthis algorithm may depend on the initial conditions and maybe considerably slower than that of AP algorithms based onpure projections [49]. D. Computational complexity
Except for the projection operator P S $ , each iteration ofthe three formulated AP algorithms relies on vector addition,scalar product, and value update. These are linear in thenumber of sample points. The operator P S $ can be easilyimplemented by using the direct and inverse fast Fourier trans-forms (FFTs) and setting the relevant elements of the signal tozero in the Fourier domain. The current state-of-the-art FFTalgorithms feature O ( n log n ) time complexity [50], which,thus, sets the overall time complexity of the AP algorithmsintroduced in this work. Our numerical experiments suggestthat the convergence speed in terms of iteration number isindependent of the signal length (see Suppl. Mat. M).IV. P ERFORMANCE T ESTS
To evaluate the AP algorithms introduced above, we com-pared their performance with the AS and LDC demodulationapproaches when applied to infer the modulator of predefinedsynthetic test signals. The LDC approach was implementedby using two state-of-the-art quadratic programming solvers:Gurobi (v8.1.1) [51] and OSQP (v0.6.0) [52]. The AS demod-ulation was achieved by using the FFT-based approach [53].In that case, we additionally low-pass filtered the obtainedmodulator estimate with P S $ to regularize it. We refer to thismodified demodulation scheme as AS-LP. A. Test signals
The test signals were composed as products of a modulatorand a carrier: s = m ◦ c . Four types of c , approximatingbasic building blocks of real-world signals, were used: non-stationary sinusoidal , harmonic , and spike-train , as well asstationary white-noise (see, respectively, (176), (180), (184),and (188) in Suppl. Mat. I). The former two were combinedwith modulators of nonstationary Gaussian origin, while thelatter two types of carriers were paired with the so-called maximally uniformly-distributed modulators (see, respectively,(160) – (162) and (163) – (168) in Suppl. Mat. I). t (s)
010 0.1 0.2 0.3 t (s)
010 0.1 0.2 0.3 t (s) t (s) x sinusoidal A white-noiseharmonic B C D spike-train
Fig. 2. Typical examples of the test signals (gray), featuring nonstationary sinusoidal ( A ), harmonic ( B ), spike-train ( C ), and stationary white-noise ( D )carriers, and their modulators obtained by using the AP-B (red) and AS-LP (green) algorithms. The signals are represented by their absolute values here. Thepredefined modulators are shown in black. In all cases, modulator and carrier pairs were selected tomeet the core recoverability condition d n/d e ≥ $ − , atleast approximately. The remaining parameters of m and c (see Suppl. Mat. I) were chosen to imitate realistic conditionsas much as possible. For example, in all cases, signals weretaken as segments of longer time series, and thus, werenot n -periodic. The center frequencies of the sinusoidal andharmonic carriers were set so that only sample points with | c i | ≈ rather than | c i | = 1 were available. B. Performance evaluation
Demodulation performance was evaluated by using twocomplementary measures: 1) error of the modulator estimate E m = k m − ˆm k / k m k ; and 2) execution time of the algo-rithm on the computer T cpu . We evaluated the AP and LDCalgorithms in the mode when T cpu depends on the total numberof sample points but not on the effective degrees of freedom.This choice made the results general, independent of a selectedcutoff frequency $ . To insure against outliers, we averaged E m and T cpu over ten independent signal realizations.Execution of the AP and LDC algorithms is controlled bya set of metaparameters whose choice influences the output.Therefore, we aimed for the Pareto fronts, not separate points,in the ( E m , T cpu ) plane. Due to the computing speed limita-tions inherent to the LDC approach, we had to exploit signaldecomposition into separate fragments for this analysis. In par-ticular, signals were split into segments that were demodulatedseparately and then put together [19]. This allowed achievinga linear growth in the computation time with the total lengthof the signal, and hence, speeding up the calculations. Afteridentifying the optimal control-parameter combinations, wecompared all methods by demodulating whole signals.Sets of the demodulation control parameters that we consid-ered for the Pareto optimality analysis, including those defin-ing the signal splitting, are provided in Suppl. Mat. J. Detailson the execution of the performance tests on a computer canbe found in Suppl. Mat. K. C. Results
Fig. 2 shows representative fragments of the test signalsfrom all four classes (gray) and their modulator estimates ob-tained by using the AS-LP (green) and AP-B (red) algorithms.Whereas the AP-B allows obtaining high-quality estimates ˆm in all four cases, the AS-LP does so only for sinusoidal signals. Results of the performance evaluation in the form of Paretofronts in the ( E m , T cpu ) plane for n = 2 are displayed inFig. 3 A–D. Panels E–H of the same figure show T cpu vs. n relations derived by using no window splitting. A closeranalysis of these data reveals the following:1) The AP algorithms feature lower bounds on the demod-ulation error E m than the LDC method (Fig. 3 A–D).2) The AP algorithms are up to five orders of magnitudefaster than their LDC counterparts for achieving thesame E m when optimal signal window splitting is used(Fig. 3 A–D). The difference is even more pronouncedwhen no window splitting is assumed (Fig. 3 E–H). Forexample, to process a 1 s length signal sampled at16 kHz, the LDC needs s of CPU time, in contrastto − s taken by the AP-A.3) T cpu varies substantially (up to three orders of magni-tude) even between different AP algorithms (Fig. 3 A–D). The AP-A ranks as the fastest, and the AP-P as theslowest one for all tested signals.4) Despite the differences in T cpu , all AP algorithms featuresimilar lower bounds on E m , except the spike-trainsignals, when the AP-B and AP-P can noticeably surpassthe AP-A on the relative scale (Fig. 3 A–D). Neverthe-less, on the absolute scale, the AP-A still performsreasonably well.5) For all tested signals, the AP-A algorithm outperformsthe AS-LP based demodulation in the sense that it canachieve the same or smaller errors with the same T cpu (Fig. 3 A–D). Moreover, compared with the AS-LP, APalgorithms exhibit much lower bounds on E m .6) Even without the window splitting, the AP-A algorithmprovides excellent demodulation quality with T cpu only2–3 times above that of the AS-LP method (Fig. 3 E–H).We found that the decrease in E m along the Pareto frontsis mainly determined by an increase in the demodulationwindow. In particular, the lower bounds on E m are achievedby the particular algorithms when the signal is demodulatedusing no window splitting. The relatively lower precision ofthe AP-A algorithm compared with AP-B and AP-P in the caseof nonstationary spike-trains can be reduced to its accelerationmechanism. Indeed, in the AP-A, upscaling of iterates a ( i ) is effectively based on the averaging of ( m ( i ) − a ( i ) ) overa window of length ≈ n/ (2 π$ ) at each sample point. Theprecision of these estimates is more vulnerable to deviationsfrom the exact recovery conditions for sparse carriers. -2 -1 E m -4 -2 -4 -3 -2 -1 E m -4 -2 -2 -1 E m -4 -2 -3 -2 -1 E m -4 -2 T c pu ( s ) AP-PAP-BAP-ALDC-GLDC-O
AS-LP n -6 -4 -2 E m = 4 10 -2 E m = 3 10 -2 E m = 4 10 -1 n -6 -4 -2 E m = 5 10 -2 E m = 2 10 -2 E m = 9 10 -1 n -6 -4 -2 E m = 2 10 -2 E m = 1 10 -2 E m = 5 10 -1 n -6 -4 -2 T c pu ( s ) LDC-GAP-A
AS-LP E m = 7 10 -3 E m = 3 10 -3 E m = 3 10 -3 sinusoidal A white-noiseharmonic B C D spike-train
GFE H
Fig. 3. Performance evaluation. A – D : Pareto fronts in the ( E m , T cpu ) plane for different demodulation algorithms applied to the four different types oftest signals when window splitting is used. Green stars mark the results of the AS-LP method. Color arrowheads point to the lower bounds on E m for therespective AP algorithms. Black arrowheads indicate demodulation error E m values of the AP-B algorithm calculated locally for signal windows shown inFig. 2. E – H : Dependence of the demodulation time T cpu on the signal length n at (cid:15) tol = 10 − when window splitting is not exploited. E m values in thelegends correspond to demodulation results at n = 2 ≈ . · (see color arrowheads). i -20 -15 -10 -5 i -20 -15 -10 -5 i -20 -15 -10 -5 i -20 -15 -10 -5 AP-BAP-PAP-A 10 i -2 -1 sinusoidal A white-noiseharmonic B C D spike-train i -3 -2 -1 i -2 -1 i -3 -2 -1 E m AP-B AP-PAP-A
GFE H
Fig. 4. Convergence analysis of the AP algorithms. A – D : Dependence of the infeasibility error (cid:15) on the iteration number i for the AP-B, AP-A, and AP-Palgorithms applied to the four different types of test signals with n = 2 and no window splitting. E – H : Analogous plots to A–D made for the demodulationerror E m instead of (cid:15) . Dotted lines show hypothetical E m values that would be obtained if we continued the AP-A iterations after reaching the final solution. As can be expected, the high accuracy of modulator esti-mates achieved by the AP algorithms implies the high qualityof carrier predictions ˆc = s / ˆm (see Suppl. Mat. L and Fig. 11therein). The AP approach leaves the AS-LP behind in termsof carrier estimation for all four signal types considered (seeSuppl. Mat. L). When applicable, the inferred ˆc can be furtherfrequency-demodulated by using dedicated techniques (see [1],[54] and references given there).The impressive performance of the AP-A algorithm in termsof E m , E c , and T cpu makes it an ideal candidate for amplitudedemodulation of arbitrary bandwidth signals. Its AP-B and AP-P counterparts can be used instead if higher precision is neededin specific cases, as illustrated by the spike-train signals above. V. C ONVERGENCE T ESTS
To clarify the differences between the T cpu estimates of thethree AP algorithms and understand the relationship betweenthe demodulation and infeasibility errors, we performed aconvergence analysis with the test signals from the previoussection. The simulation results for fixed n = 2 using nowindow splitting are summarized in Fig. 4. A closer inspectionuncovers the following:1) The convergence rates in terms of both (cid:15) and E m parallelthe differences in the computing speed of different APalgorithms. Among them, the fastest is the AP-A, whichreaches any given (cid:15) or E m level with the smallest num-ber of iterations. The AP-P algorithm is the slowest one. P(0) = 0.25 t (s) P(0) = 0.65 t (s) P(0) = 0.90 t (s) x P(0) = 0.70 t (s)
01 0 0.25 0.5 0.75 1
P(0) -2 -1 P(0) -3 -2 -1 P(0) -3 -2 -1 P(0) -3 -2 -1 E m AP-BAP-PAP-AAS-LP sinusoidal A white-noiseharmonic B C D spike-train
GFE H
Fig. 5. Robustness evaluation. A – D : Dependence of the demodulation error E m on P (0) (the probability of missing points) for the four types of test signalsand different AP algorithms at (cid:15) tol = 10 − (color coding). E – H : Representative examples of demodulation at various P (0) levels for the test signals fromA–D. Color code: gray – the absolute-value signal, black – the original modulator, color – modulators inferred by different algorithms.
2) The AP-A algorithm converges in a finite number ofiterations ( < ) for all types of test signals studied. Inparticular, it requires only ≤ iterations to reach theplateau level of the demodulation error E m . This factexplains the extraordinary computational efficiency ofthe AP-A documented in Section IV.3) Differently from the convergence error k m ( i ) − m † k / √ n , the dependence of E ( i ) m on i can be non-monotonic if m † is not strictly equal to m (Fig. 4 E, G).Then, E ( i ) m starts growing with increased i after reachingthe minimum point. However, this growth is mild andof no practical importance as long as m † ≈ m , i.e., ˆm ≈ m .The results shown in Fig. 4 represent only signals of fixedlength ( n = 2 sample points). Additional simulations sug-gested no dependence on n (see Suppl. Mat. M).VI. R OBUSTNESS T ESTS
The pivotal condition for successfully separating themodulator-carrier information of a given signal by our ap-proach is d n/d e ≥ ω − . In practice, this requirementis not necessarily met. Hence, the choice of a particulardemodulation method must be guided not only by the algo-rithmic efficiency but also robustness to deviations from theideal recovery conditions. To shed light on this aspect, weconsidered demodulation of the test signals from Section IV-Acorrupted by a multiplicative Bernoulli- { , } noise. In thissetup, sample points, including the decisive | s i | = m i , areeliminated with the probability of “0” elements of the noise( P (0) ), effectively decreasing the value of d n/d e .We found that all three AP algorithms considered in thiswork show a similar degree of robustness to increased P (0) (see Fig. 5). Only in the case of sinusoidal signals, the AP-Ais slightly inferior to the AP-B and AP-P. Interestingly, theadvantage of the AP-B and AP-P over the AP-A in the caseof spike-train signals discussed in Section IV-C disappearsin the presence of even small distortions (see Fig. 5 C). The differences in the E m vs. P (0) relations seen in Fig. 5 A–D are predetermined by different densities of | c i | ’ pointsinherent to each carrier type. Analogous results to those shownin Fig. 5 A–D are obtained when considering carrier recoveryvia ˆc = s / ˆm (see Fig. 12 in Suppl. Mat. L).In contrast to the AP approach, the AS-based demodulationis highly vulnerable to missing sample points, and hence, todecreased d n/d e (Fig. 5 A–D). Even for sinusoidal signals,which the AS and AS-LP are specially designed for, thezeroing of data points leads to a rapid decline in demodulationquality (Fig. 5 A, E).The robustness to missing sample points endows the AP de-modulation method with a highly valuable practical advantage.In particular, it can be exploited in real-world situations when:1) the sampling rate is low; 2) some segments of the signalvalues are lost; 3) some sample points are corrupted by noisesuch that the level of these points can be reduced below the realmodulator by low-pass filtering or explicitly identifying them.In this context, the PAD and LDC demodulations compare tothe AP approach by construction [23].VII. H IGH -L EVEL P ROPERTIES
As emphasized in Section II, different demodulation meth-ods can be derived by requesting adherence of the inferredmodulators and carriers to a set of particular properties.Typically, various combinations that consist of a few out ofmany reasonable requirements are sufficient for unique demod-ulation formulations. However, some of these requirements areinconsistent with each other, making virtually all classical de-modulation approaches fail to satisfy one or another essentialcondition [9], [23], [25]. For example, the AS demodulationmethod may return an unbounded modulator estimate for abounded signal [25].The AP approach formulated in this work is compatiblewith the following high-level requirements, which have crys-tallized as inseparable from the notion of proper amplitudedemodulation with time [19], [23, Section 3.5.2]: = 111 Hz = 222 Hz f c = 800 Hz t (s) x = 22.1 kHz = 222 Hz= 55 Hz t (s) x "…at one hundred hertz…" AB Fig. 6. Direct demodulation of speech signals. A : A band-pass-filtered signal of an utterance “. . . at one hundred hertz . . . ” by a female speaker; the signalwas obtained with an equivalent rectangular bandwidth filter of the cochlea centered at f c = 800 Hz and ∆ = 111 Hz . B : The original signal of the utteranceused in panel A (full bandwidth of . ). In both panels, gray color marks the absolute-value version of the signals considered for demodulation. Blueand red lines show their modulators obtained by using the AP-A algorithm with the cutoff frequency ω set to, respectively,
55 Hz and
222 Hz . Insets displaytime-expanded segments of the original window. Red arrowheads indicate the ringing artifacts of the modulator at some prolonged intervals of low signallevels. Source of the original signal: audio edition of
The Economist magazine, issue March 19th 2016, article “Restoring lost memories. • Boundedness:
The modulator and carrier of a boundedsignal are bounded. In particular, it is required that −∞ The modulator and carrier of a scaledsignal are equal to the modulator and carrier obtainedfrom the original signal and then scaled by the sameamount. The adherence of the AP approach to this con-dition follows from two facts. First, projection operators P S ≥| s | and P S $ are homogeneous with degree 1, i.e., P S [ α s ] = α P S [ s ] . Second, each iteration of the APalgorithms can be expressed as a weighted sum of theseprojections with the weights independent of the scale. • Smoothness: The modulator of a bounded signal in itscontinuous-time representation is smooth. Because weuse a discrete-time representation, this requirement has tobe adjusted. In particular, let us denote by m t and m t +∆ t the finite-difference approximations of the modulator’stime-derivatives of any order at two subsequent timepoints t < t + ∆ t . Then, we require that, for any (cid:15) > ,there exists a δ > such that | m t +∆ t − m t | < (cid:15) when | ∆ t | < δ . The AP approach satisfies this requirementthrough the boundedness of the modulator and the band-width constraint set by S $ on it. • Idempotence: Information associated with the qualitiesof modulators and carriers is fully separated. Specifi-cally, demodulation reapplied to an estimated modulator(carrier) must return the same modulator (carrier). TheAP approach satisfies the idempotence requirement forthe modulator exactly. Indeed, when any AP algorithm is applied to its final solution ˆm = m † , the latteris recognized as the final solution again after the firstnew iteration by construction. Regarding the carrier, theidempotence holds whenever the recovery conditions dis-cussed in Section II-C are met. That is because, in thosecases, ˆc resulting from the first demodulation contains asufficient number of | c i | = 1 points to uniquely define the m = as the norm-minimizing element of S ≥| ˆc | ∩ S $ .If the recovery conditions are met only approximately,we expect no marked deviations from the idempotencecondition (see Fig. 14 in Suppl. Mat. N).By fulfilling the above requirements, the AP approach parallelsthe methods of PAD and LDC demodulation [19], [23]. In thissense, all of them outperform the classical techniques.VIII. D EMODULATION OF S PEECH S IGNALS Amplitude demodulation is of central importance in var-ious tasks of processing and analysis of speech signals.Application-wise, this procedure is used in hearing restoration[10], [55], speech recognition [2], [56], [57], and source sepa-ration [58], [59]. On the theory side, amplitude demodulationis exploited in neurophysiological and psychophysical studiesof auditory information processing in the brain [11], [12],[60], [61]. Depending on the problem, demodulation of eithernarrow subband [56], [58], intermediate subband [10], [11], orwhole wideband signal [12], [62] is needed. In all these cases,the modulators and carriers convey the information aboutspecific aspects of speech, e.g., semantic meaning, associatedemotion, or speaker identity, that need to be extracted.In this section, we apply the newly-introduced AP approachto speech demodulation to further demonstrate its potential.To represent the range of possible real-world situations, weconsider two limiting signal types: 1) a narrow subband = 16 kHz= 55 Hz t (s) x "…protein which forms p…" Fig. 7. Demodulation of speech signals using dynamic range compression. An audio signal of “. . . protein which forms p . . . ” uttered by a female speaker(full bandwidth of 16 kHz ). Gray – the absolute-value version of the signal considered for demodulation, red – its modulator obtained by using the AP-Aalgorithm with ω = 55 Hz and no compression (as in Fig. 6), black – modulator of the signal obtained when employing the dynamic range compression [see(13)], violet – interpolation of the latter two [see (14)]. Red arrowheads indicate ringing artifacts of the modulator estimate. Source of the original signal: thesame as Fig. 6. component of a signal obtained by a standard auditory ERBfilter [63] and 2) an original wideband signal. A. Direct demodulation By construction, the output of auditory ERB filters occupiesa frequency subband whose width ∆ is much smaller than itscenter frequency f c [63]. The resulting signal is an amplitude-and phase-modulated sinusoidal s = m ◦ sin(2 πf c t + ϕ ) , withmost of the energies of m and ϕ residing in the frequencyinterval [0 , ∆] [64]. Hence, by the recovery conditions of theAP approach (see Section II-C), setting the cutoff frequency $ between ∆ and f c necessarily results in accurate estimates ˆm and ˆc . In particular, note that the local maximums of | s | correspond to the true sample points | s i | = m i . Thus, the highquality of demodulation is visually conveyed by a tight matchof ˆm and | s | at these sample points. We illustrate our claims inFig. 6 A, where a band-pass component of a female utterance “. . . at one hundred hertz . . . ” with f c = 800 Hz , ∆ = 111 Hz ,and $ = 222 Hz is considered. Taking into account thatthe local maximum points | s i | are locally regular and thatthey correspond to m i , we could exploit Proposition A.2 inSuppl. Mat. A to find that E m ≤ · − .Wideband speech signals are more challenging than theirnarrow subbands. They are built of temporarily structured seg-ments of quasi-random and quasi-harmonic carriers, possiblyfeaturing frequency glides [65]. These carriers are amplitude-modulated at different timescales, ranging between a hun-dred milliseconds and several seconds [23], [66]. The powerspectral density of the corresponding modulators is vanish-ingly small above 20 Hz (see Fig. 1 in [67]). Moreover, aswe demonstrate in Suppl. Mat. O, the carrier components ofnatural speech signals align to the recoverability conditions ofthe AP approach for m ∈ M ω with ω up to at least ∼ 50 Hz .Therefore, we expect appropriate performance from the APalgorithms in the setting of wideband speech.Fig. 6 B displays demodulation results of the full-band ver-sion of the speech segment considered in Fig. 6 A by the AP-Aalgorithm with $ = 55 Hz . The obtained ˆm (red) envelopsseparate phonemes of the sound waveform tightly, indicatingappropriate recovery of the true m (see Section VIII-B next).However, intervals corresponding to prolonged transitions be-tween phonemes or words are corrupted by ringing artifacts (marked by red arrowheads in Fig. 6 B), implying the necessityof higher frequency components to represent these transitions.Hence, although the power spectral density of the true m isvery low above 20 Hz , it sums to a noticeable contribution.Unfortunately, any attempt to cancel the artifacts by justincreasing $ fails by breaking the recovery conditions, asillustrated by the blue line in Fig. 6 B ( $ = 222 Hz there).No improvement is achieved by utilizing the AP-B, AP-P, orLDC algorithms either (data not shown). B. Demodulation using dynamic range compression The aforementioned problem with modulator estimates ofsignals with sharp transitions to/from prolonged intervals oflow-signal amplitude can be resolved by using a dynamicrange compression. In particular, instead of demodulating theoriginal signal s directly, we first apply a chosen AP algorithmto its compressed version: s = sgn( s ) ◦ | s | /p . (12)Here, p ∈ (1 , + ∞ ) controls the level of compression. Themodulator estimate ˆm ∗ of s is then evaluated by inverse-transforming the modulator ˆm of s : ˆm ∗ = ˆm p . (13)The idea behind (12) is that the compression makes signalsmore uniform and, effectively, smooths their sharp changesresponsible for ringing artifacts in the modulator estimates.These sharp changes are restored in the modulators withoutartifacts by the inverse transform (13).The expected effect of the compression procedure is illus-trated in Fig. 7, where signal demodulation of an utterance “. . . protein which forms p. . . ” is considered. Differently fromthe direct demodulation result ˆm (red line), the estimate ˆm ∗ obtained by using the compression with p = 3 (black line)shows good alignment with | s | in the segments of both lowand high intensity. To justify that this alignment really reflectsthe recovery of the true modulator, we performed additionaltests where chimeric signals built of ˆm ∗ from Fig. 7 andnatural speech carriers were demodulated (see Suppl. Mat. O).We found low demodulation errors, with E m ranging between · − and · − for different carrier components of speechsignals (see Fig. 16). = 111 Hz f c = 800 Hz = 333 Hz T cpu = 1.6 %= 6 ms t (s) x = 8 kHz= 40 Hz T cpu = 3.2 %= 48 ms t (s) x A "…with little human hand-holding…" B Fig. 8. Demodulation in real-time. A : A band-pass-filtered signal of an utterance “. . . with little human hand-holding . . . ” by a male speaker; the signal wasobtained with an equivalent rectangular bandwidth filter of the cochlea centered at 800 Hz (bandwidth of 111 Hz ). B : The original signal of the utteranceused in panel A (full bandwidth of ). In both panels, gray color marks the absolute-value version of the signals considered for demodulation. Blue andred lines show their modulators obtained by using the real-time version of the AP-A algorithm with the cutoff frequency ω set to, respectively, 333 Hz and 40 Hz . Source of the original signal: audio edition of The Economist magazine, issue March 19th 2016, article “Artificial intelligence and Go.” The compression level p = 3 used above was adjusted bya trial and error for speech signals. In general, the gains inaccuracy at low levels with increased p comes at the expenseof reduced precision of modulator estimated at high signallevels. Thus, a compromise between those two effects must bereached to find an optimal p . Moreover, the precision of themodulator estimates can be further increased by interpolatingbetween ˆm ∗ (more accurate for low signal levels) and ˆm (more accurate for high signal levels). For example, the violetline in Fig. 7 shows a weighted average of the form ˆm (cid:5) = ˆm ◦ w + ˆm ∗ ◦ (1 − w ) , (14)where w i = (cid:18) − e a · ( ˆm ∗ i / max[ ˆm ∗ ]) e a · ( ˆm ∗ i / max[ ˆm ∗ ]) − b (cid:19) · (cid:18) − e a e a − b (cid:19) − (15)for i ∈ I n , with b = 3 and a = 10 . In general, anoptimal interpolation between ˆm and ˆm ∗ can be learned byminimizing k ˆm (cid:5) k over a chosen class of functions. Othercompression models than (12), e.g., s = sgn( s ) ◦ log(1+ p ·| s | ) ,can be used to evaluate ˆm ∗ as well. C. Demodulation in real-time A number of amplitude demodulation applications, e.g.,speech recognition [2], ultrasound imaging [68], and cochlearprosthesis [55], necessitate real-time processing. As wedemonstrate below, the exceptional computational efficiencyof the AP approach allows it to fulfill that requirement.The nature of the task implies that online modulator es-timates have to be generated by sequentially demodulatingwindowed segments s ( j ) of a signal s at each updated samplepoint j across time: s ( j ) : s ( j ) i = w i · s j − k l − i i ∈ { , , . . . , k } . (16) Here, k is the number of sample points corresponding to thesegment, and k l denotes the number of sample points of it thatare to the left of the current point j . w i , w , . . . , w k are vectorelements of the window function. The real-time modulatorestimate ˆ m ?j at sample point j is calculated as ˆ m ?j = ˆ m ( j ) k l +1 , (17)where ˆm ( j ) is a modulator estimate of s ( j ) .It follows from the time-frequency uncertainty principle [8]that accurate evaluation of ˆ m ?j requires s ( j ) with a durationof the order of the inverse of the effective bandwidth of themodulator, or longer. This condition sets the lower bounds onthe segment length k and sampling delay k τ = k − k l − of ˆm ? . We found empirically that k τ ≈ · ( f s /$ ) and k ≈ · ( f s /$ ) are typically sufficient for accurate demod-ulation of wideband speech. These numbers are around twotimes smaller for narrow frequency band components of thesesignals. We know that $ ≥ 40 Hz for the wideband speechand its subbands. Thus, delays k τ ≤ 50 ms for estimating ˆm ? are sufficient without a sacrifice in precision then. The mainrequirement for the window function in (16) is that it smoothlyscales the signal to 0 at the boundaries, with no effect at themidst. We used a modified version of the Hann window forthis purpose: w i = sin (cid:16) π · ( i − · k l (cid:17) , ≤ i ≤ k l , i = k l + 1cos (cid:16) π · ( i − k + k τ )2 · k τ (cid:17) , k − k τ + 1 ≤ i ≤ k . (18)Fig. 8 shows simulation results of real-time demodulationof a male utterance “. . . with little human hand-holding . . . ” (sampling rate f s = 16 kHz ) based on the AP-A algorithm.There, demodulation was performed with k = 1536 and k τ = 768 ( τ = 48 ms ) for the original signal (Fig. 8 B). 100 1 10 t (s) t (s) x A B m o d u l a t i o n d e m o d u l a t i o n original × inferred × Fig. 9. Demodulation in 2D by using the AP-A algorithm. A : Synthetic fringe pattern of moderate bandwidth. Left – the pattern to demodulate, upper right –the original modulator and carrier, lower right – the inferred modulator and carrier. B : Wideband signal consisting of randomly-placed spikes of finite width.Black grid – the signal ¯s to demodulate, white grid – the original modulator m of the signal, color surface – the estimated modulator. Its subband component centered at 800 Hz (Fig. 8 A) wasprocessed with k = 128 and k τ = 65 ( τ = 4 ms ). Ineach case, ˆ m ?j was updated with a frequency of · $ . Theobtained estimates ˆm ? are in very good agreement with ˆm ∗ derived by using offline demodulation of the whole signal,with k ˆm ? − ˆm (cid:5) k / k ˆm (cid:5) k < . . Importantly, they wereachieved with modest CPU usage: T cpu amounted to only1.6 % (subband signal) and 3.2 % (wideband signal) of thetime length of the demodulated signal on an Intel Core i7-7700 CPU run in single-thread mode. For comparison, thesenumbers were, respectively, ∼ · and ∼ · timeshigher for the LDC method.An advantageous side effect of splitting the signal into smallwindows for demodulation is that it prevents the ringing arti-facts (compare Fig. 8 B and Fig. 6 B). This is so because signallevels do not typically spread over different scales in a shorttime window. The window splitting also allows generalizingdemodulation to situations when the cutoff frequency ω of themodulator varies broadly.IX. E XTENSIONS AND G ENERALIZATIONS A. Demodulation in higher dimensions Amplitude demodulation has found successful applicationsbeyond the setting of 1D signals. Several 2D extensions ofthe classical AS approach have been introduced and usedfor solving tasks in computer vision [7], [69], analysis ofspeech spectrograms [70], and biomedical imaging [4], [71],[72]. These methods are limited to locally narrowband signals,which manifest visually as fringe patterns (see Fig. 9 A). Thisbandwidth restriction is evaded by a generalization of theAP approach to higher dimensions that we present next. Theextension is immediate and follows from intuitive abstractionsof the constraint sets introduced in Section II.Consider a D -dimensional signal s ( t , t , . . . , t D ) . Itsuniformly sampled version s is an element of an n -dimensional Euclidean space T nD of real-valued order D tensors with n = Q Di =1 n i and an inner product h s (1) , s (2) i = P n i =1 · · · P n D i D =1 (¯ s (1) i ··· i D · ¯ s (2) i ··· i D ) . The respec-tive D -dimensional DFT is given by F = F (1) ⊗ F (2) ⊗ · · · ⊗ F ( D ) , (19) where F ( i ) is a unitary DFT defined over R n i . Then, theanalogs of the constraint sets S ≥ , S ω , S ≥| s | , S | .. |≤ , and S { } ,d from Section II read as S ≥ ¯0 = { x ∈ T nD : x i ··· i D ≥ , i j ∈ I n j } , S ω = { x ∈ T nD : ( Fx ) i ··· i D = 0 , i j ∈ ( I n j \ I ω j n j ) } , (20) S ≥| ¯s | = { x ∈ T nD : x i ··· i D ≥ | ¯ s i ··· i D | , i j ∈ I n j } , and S | .. |≤ ¯1 = { x ∈ T nD : | x i ··· i D | ≤ , i j ∈ I n j } , S { } , d = (cid:8) x ∈ T nD : ( ∀ i ··· i D ) R ( i ··· i D , x , d ) ≥ , ( ∃ i ··· i D ) R ( i ··· i D , x , d ) = 1 } , (21)where R ( i ··· i D , x , d ) = P j ≥ i ··· P j D ≥ i D (cid:2) I { } ( | x j ··· j D | ) · θ (cid:0) − P Dk =1 ( i k − j k ) /d k (cid:1)(cid:3) − I { } ( | x i ··· i D | ) . (22)Simply substituting (20) – (21) for their D = 1 versions in(2), (4), (6), (8), (9), and (10) generalizes the modulator M ω and carrier C d sets as well as the modulator estimator ˆm andthe respective AP algorithms. In particular, an m ∈ M ω is anonnegative signal with a low-pass rectangular spectrum setby ω = ( ω , . . . , ω D ) along each of the D dimensions in theDFT domain. An c ∈ C d is a signal bounded between − and with the | ¯ c i ··· i D | = 1 sample points packed sufficientlydensely, as implied by d = ( d , . . . , d D ) .Without providing formal proofs, we state that all proposi-tions and assertions of Sections II and III about the modulatorrecoverability and convergence of the AP algorithms general-ize to D -dimensional signals defined above. All quantitativeconditions involving the parameters $ , ω , d , n , and n s in the D = 1 case are then replaced by elementwise conditions for $ i , ω i , d i , n i , and n s,i at i ∈ I D .Fig. 9 illustrates the potential of the AP-A algorithm withthe help of two D = 2 cases. Fig. 9 A shows successfuldemodulation results for a synthetic narrowband fringe pattern( E m = 1 · − , E c = 3 · − ). Fig. 9 B displays high-accuracydemodulation of a wideband signal built of randomly-placedspikes of finite width as c and a Gaussian random field witha rectangular amplitude spectrum as m . There, the white grid corresponds to the original m , while the color surfacerepresents its estimate ˆm ( E m = 4 · − , E c = 1 · − ).The ability of the AP approach to deal with widebandsignals allows it to cover a wider range of practically relevantsituations. Among examples are nonlinear ultrasound imaging[4], [14], speech processing [20], [70], and complicated casesof optical interference/diffraction setups [73]. Moreover, itcan also be of great use in time-critical imaging settings byproviding high modulator estimation accuracy at low samplingrates of the signal (see, e.g., [72]).The minimum number of sample points necessary to coversimultaneously for appropriate demodulation increases expo-nentially with D . Therefore, the computational advantage ofthe AP over the PAD and LDC demodulation approachesis even more pronounced in higher dimensions. In fact, ifevaluated by using the FFT method, F features a O (cid:0) n log n (cid:1) computational time complexity. Hence, the time complexity ofthe AP algorithms is defined by the total number of samplepoints of the signal irrespective of its dimensionality. B. Generalized modulators and nonuniform sampling The demodulation approach formulated in the present workbuilds on the assumption that modulators are nonnegative ele-ments of a low-pass DFT subspace of T nD . However, as followsfrom the convergence proofs in Suppl. Mat. F, all introducedAP algorithms are bound to converge to an ˆm ∈ M ω and ˆc ∈ C d independent of the origin of the linear subspacebehind M ω . This naturally raises a question of whether theAP algorithms could recover true m and c under a generalizedsubspace assumption. Our preliminary experiments suggesta positive answer but subject to extra recovery conditionsspecific to a subspace of choice.For example, consider a subset of ω − randomly chosenbasis vectors of the DFT over R n . Denote the correspondingspace as F ω . It can be shown by example that a systemresulting from random subsampling of aforementioned vectorsat ω − time points may be linearly dependent. If so, it thenfollows from the proof of Proposition II.1 that, in contrast toan m ∈ S ω , full recovery of an m ∈ F ω necessitates morethan ω − true sample points.The problem of formulating modulator recovery conditionsfor different linear subspaces sets directions for future studies.If successful, these extensions would allow to:1) broaden the concept of the amplitude modulator beyondthe low-pass DFT signals,2) loosen the constraints on the positioning of the | c i | = 1 sample points for recoverable carriers whenever a morecompact representation of modulators is available,3) encompass nonuniform sampling.While the above points are yet to be demonstrated rigor-ously, the results of the present work already provide a strategyfor an arbitrarily-approximate nonuniform sampling. Indeed,for any time grid ˜t ∈ R ˜ n , we can find a uniform grid t ∈ R n such that, for every j ∈ I ˜ n , there exists an i ∈ I n with | ˜ t j − t i | being arbitrarily small. We can then interpolate theoriginal data ˜s ∈ R ˜ n on the uniform grid t by s i = ( ˜ s j , | ˜ t j − t i | = min[ | ˜ t j · − t | ]0 , otherwise , i ∈ I n , (23)to obtain an s ∈ R n . The bandwidth constraint on m impliesthat all components of s corresponding to the true samplepoints of ˜s are desirably close to the true sample points of theoriginal signal if n is large enough. Then, Proposition II.2 as-sures that modulator-carrier recovery is possible via (6) underthe conditions discussed in Section II-C for uniformly sampledsignals. The described strategy requires increasing the effectivedimensionality of the signal. This may still be more efficientthan evaluating metric projections onto subspaces spanned byarbitrary nonuniform sampling basis vectors, which are notorthogonal in general.X. C ONCLUSION In this paper, we introduce a new approach to amplitudedemodulation of arbitrary-bandwidth signals. We frame de-modulation as a problem of signal recovery from an unlabeledmix of its true and corrupted sample points. A judicious choiceof structural assumptions on the modulator and carrier enableshigh-quality separation of the latter two for a wide range ofrelevant cases. We show that this goal can be reached through aconstrained norm minimization of the modulator and formulatetailor-made alternating projection algorithms to achieve that inpractice.The generality and numerical efficiency of the new approachmake it a preferred choice in many situations. In the context ofnarrowband signals, the new method outperforms the classicalalgorithms in terms of robustness to data distortions andcompatibility with nonuniform sampling. When consideringthe demodulation of wideband signals, it surpasses the currentstate-of-the-art techniques in terms of computational efficiencyby up to many orders of magnitude. Such performance enablespractical applications of amplitude demodulation in previouslyinaccessible settings. Specifically, online and large-scale of-fline demodulation of wideband signals, signals in higherdimensions, and poorly-sampled signals become practicallyfeasible. The algorithms underlying the new approach aresimple and easy to code.A CKNOWLEDGMENT The author thanks his colleagues K. Husz´ar and G. Tkaˇcikfor valuable discussions and comments on the manuscript.R EFERENCES[1] D. Vakman, Signals, oscillations, and waves: A modern approach .Artech House, 1998.[2] B. E. D. Kingsbury, N. Morgan, and S. Greenberg, “Robust speechrecognition using the modulation spectrogram,” Speech Commun. ,vol. 25, no. 1, pp. 117–132, 1998.[3] M. G. Ruppert, D. M. Harcombe, M. R. P. Ragazzon, S. O. R.Moheimani, and A. J. Fleming, “A review of demodulation techniquesfor amplitude-modulation atomic force microscopy,” Beilstein J. Nan-otechnol. , vol. 8, no. 1, pp. 1407–1426, 2017.[4] C. Wachinger, T. Klein, and N. Navab, “The 2D analytic signal for en-velope detection and feature extraction on ultrasound images,” MedicalImage Analysis , vol. 16, no. 6, pp. 1073–1084, 2012. [5] P. Y. Ktonas and N. Papp, “Instantaneous envelope and phase extractionfrom real signals: Theory, implementation, and an application to EEGanalysis,” Elsevier Signal Process. , vol. 2, no. 4, pp. 373–385, 1980.[6] M. T. Taner, F. Koehler, and R. E. Sheriff, “Complex seismic traceanalysis,” Geophysics , vol. 44, no. 6, pp. 1041–1063, 1979.[7] K. G. Larkin, D. J. Bone, and M. A. Oldfield, “Natural demodulation oftwo-dimensional fringe patterns. I.” J. Opt. Soc. Am. A , vol. 18, no. 8,pp. 1862–1870, 2001.[8] D. Gabor, “Theory of communication. part 1: The analysis of informa-tion,” J. Inst. Elec. Eng. Part III , vol. 93, no. 26, pp. 429–441, 1946.[9] D. Vakman, “On the analytic signal, the Teager-Kaiser energy algorithm,and other methods for defining amplitude and frequency,” IEEE Trans.Signal Process. , vol. 44, no. 4, pp. 791–797, 1996.[10] B. S. Wilson, C. C. Finley, D. T. Lawson, R. D. Wolford, D. K.Eddington, and W. M. Rabinowitz, “Better speech recognition withcochlear implants,” Nature , vol. 352, no. 6332, pp. 236–238, 1991.[11] Z. M. Smith, B. Delgutte, and A. J. Oxenham, “Chimaeric sounds revealdichotomies in auditory perception,” Nature , vol. 416, no. 6876, pp. 87–90, 2002.[12] U. Goswami, “Speech rhythm and language acquisition: An amplitudemodulation phase hierarchy perspective,” Ann. N. Y. Acad. Sci. , vol.1453, no. 1, pp. 67–78, 2019.[13] S. Lin, “Demodulating wide-band ultrasound signals,” U.S. PatentUS6 248 071B1, 2001.[14] F. A. Duck, “Nonlinear acoustics in diagnostic ultrasound,” UltrasoundMed. Biol. , vol. 28, no. 1, pp. 1–18, 2002.[15] G. L. Gottlieb and G. C. Agarwal, “Filtering of electromyographicsignals,” Am. J. Phys. Med. Rehab. , vol. 49, no. 2, p. 142, 1970.[16] J. Felblinger and C. Boesch, “Amplitude demodulation of the electro-cardiogram signal (ECG) for respiration monitoring and compensationduring MR examinations,” Magn. Reson. Med. , vol. 38, no. 1, pp. 129–136, 1997.[17] D. Gill, N. Gavrieli, and N. Intrator, “Detection and identificationof heart sounds using homomorphic envelogram and self-organizingprobabilistic model,” in Computers in Cardiology , 2005, pp. 957–960.[18] W. Liu and B. Santhanam, “Wideband image demodulation via bi-dimensional multirate frequency transformations,” J. Opt. Soc. Am. A ,vol. 33, no. 9, pp. 1668–1678, 2016.[19] G. Sell and M. Slaney, “Solving demodulation as an optimizationproblem,” IEEE Audio, Speech, Language Process. , vol. 18, no. 8, pp.2051–2066, 2010.[20] ——, “The information content of demodulated speech,” in IEEE Proc.ICASSP’10 , 2010, pp. 5470–5473.[21] R. Libbey, Signal & image processing sourcebook . Springer, 1994.[22] R. S. Platt, E. A. Hajduk, M. Hulliger, and P. A. Easton, “A modifiedbessel filter for amplitude demodulation of respiratory electromyo-grams,” J. Appl. Physiol. , vol. 84, no. 1, pp. 378–388, 1998.[23] R. E. Turner, “Statistical models for natural sounds,” Ph.D. dissertation,University College London, 2010. [Online]. Available: http://discovery.ucl.ac.uk/19231/[24] R. E. Turner and M. Sahani, “Demodulation as probabilistic inference,” IEEE Audio, Speech, Language Process. , vol. 19, no. 8, pp. 2398–2411,2011.[25] P. J. Loughlin and B. Tacer, “On the amplitude- and frequency-modulation decomposition of signals,” J. Acoust. Soc. Am. , vol. 100,no. 3, pp. 1594–1601, 1996.[26] L. Cohen, P. Loughlin, and D. Vakman, “On an ambiguity in thedefinition of the amplitude and phase of a signal,” Elsevier SignalProcess. , vol. 79, no. 3, pp. 301–307, 1999.[27] J. von Neumann, Functional operators. Vol. II: The geometry of or-thogonal spaces , ser. Annals of Mathematics Studies 22. PrincetonUniversity Press, 1951.[28] R. Escalante and M. Raydan, Alternating projection methods . SIAM,2011.[29] H. H. Bauschke and J. M. Borwein, “On projection algorithms forsolving convex feasibility problems,” SIAM Rev. , vol. 38, no. 3, pp.367–426, 1996.[30] A. Ahmed, B. Recht, and J. Romberg, “Blind deconvolution usingconvex programming,” IEEE Trans. Inf. Theory , vol. 60, no. 3, pp. 1711–1732, 2014.[31] S. Ling and T. Strohmer, “Self-calibration and biconvex compressivesensing,” Inverse Probl. , vol. 31, no. 11, p. 115002, 2015.[32] Y. Chi, “Guaranteed blind sparse spikes deconvolution via lifting andconvex optimization,” IEEE J. Sel. Topics Signal Process. , vol. 10, no. 4,pp. 782–794, 2016. [33] Y. Xie, M. B. Wakin, and G. Tang, “Simultaneous sparse recovery andblind demodulation,” IEEE Trans. Signal Process. , vol. 67, no. 19, pp.5184–5199, 2019.[34] A. Oppenheim and J. Lim, “The importance of phase in signals,” Proc.IEEE , vol. 69, no. 5, pp. 529–541, 1981.[35] H. Trussell and M. Civanlar, “The feasible solution in signal restoration,” IEEE Trans. Acoust., Speech, Signal Process. , vol. 32, no. 2, pp. 201–212, 1984.[36] D. Kundur and D. Hatzinakos, “A novel blind deconvolution scheme forimage restoration using recursive filtering,” IEEE Trans. Signal Process. ,vol. 46, no. 2, pp. 375–390, 1998.[37] Y. Yang, N. P. Galatsanos, and H. Stark, “Projection-based blinddeconvolution,” J. Opt. Soc. Am. A , vol. 11, no. 9, pp. 2401–2409, 1994.[38] P. J. S. G. Ferreira, “Iterative and noniterative recovery of missingsamples for 1-D band-limited signals,” in Nonuniform sampling: Theoryand practice , F. Marvasti, Ed. Springer, 2001, pp. 235–281.[39] L. G. Gubin, B. T. Polyak, and E. V. Raik, “The method of projectionsfor finding the common point of convex sets,” USSR Comput. Math. &Math. Phys. , vol. 7, no. 6, pp. 1–24, 1967.[40] D. C. Youla and H. Webb, “Image restoration by the method of convexprojections: Part 1– theory,” IEEE Trans. Med. Imag. , vol. 1, no. 2, pp.81–94, 1982.[41] C. Franchetti and W. Light, “On the von Neumann alternating algorithmin Hilbert space,” J. Math. Anal. Appl. , vol. 114, no. 2, pp. 305–314,1986.[42] W. B. Gearhart and M. Koshy, “Acceleration schemes for the methodof alternating projections,” J. Comput. Appl. Math. , vol. 26, no. 3, pp.235–249, 1989.[43] H. H. Bauschke, F. Deutsch, H. Hundal, and S.-H. Park, “Acceleratingthe convergence of the method of alternating projections,” Trans. Amer.Math. Soc. , vol. 355, no. 9, pp. 3433–3461, 2003.[44] R. L. Dykstra, “An algorithm for restricted least squares regression,” J.Amer. Statist. Assoc. , vol. 78, no. 384, pp. 837–842, 1983.[45] J. P. Boyle and R. L. Dykstra, “A method for finding projectionsonto the intersection of convex sets in Hilbert spaces,” in Advances inOrder Restricted Statistical Inference , ser. Lecture Notes in Statistics.Springer, 1986, pp. 28–47.[46] E. Birgin and M. Raydan, “Robust stopping criteria for Dykstra’salgorithm,” SIAM J. Sci. Comput. , vol. 26, no. 4, pp. 1405–1414, 2005.[47] F. Deutsch and H. Hundal, “The rate of convergence of Dykstra’scyclic projections algorithm: The polyhedral case,” Numer. Funct. Anal.Optim. , vol. 15, no. 5-6, pp. 537–565, 1994.[48] F. Deutsch, “Dykstra’s cyclic projections algorithm: The rate of conver-gence,” in Approximation Theory, Wavelets and Applications , ser. NATOScience Series. Springer, 1995, pp. 87–94.[49] H. H. Bauschke and J. M. Borwein, “Dykstra’s alternating projectionalgorithm for two sets,” J. Approx. Theory , vol. 79, no. 3, pp. 418–443,1994.[50] P. Duhamel and M. Vetterli, “Fast Fourier transforms: A tutorial reviewand a state of the art,” Elsevier Signal Process. , vol. 19, no. 4, pp.259–299, 1990.[51] Gurobi Optimization, LLC, Gurobi optimizer reference manual Math. Prog. Comp. ,vol. 12, no. 4, pp. 637–672, 2020.[53] L. Marple, “Computing the discrete-time “analytic” signal via FFT,” IEEE Trans. Signal Process. , vol. 47, no. 9, pp. 2600–2603, 1999.[54] R. Wiley, H. Schwarzlander, and D. Weiner, “Demodulation Procedurefor Very Wide-Band FM,” IEEE Trans. Commun. , vol. 25, no. 3, pp.318–327, 1977.[55] B. S. Wilson and M. F. Dorman, “Cochlear implants: A remarkable pastand a brilliant future,” Hear. Res. , vol. 242, no. 1, pp. 3–21, 2008.[56] S. Wu, T. H. Falk, and W.-Y. Chan, “Automatic speech emotion recog-nition using modulation spectral features,” Speech Commun. , vol. 53,no. 5, pp. 768–785, 2011.[57] B. Lee and K.-H. Cho, “Brain-inspired speech segmentation for au-tomatic speech recognition using the speech envelope as a temporalreference,” Sci. Rep. , vol. 6, p. 37647, 2016.[58] G. Hu and D. Wang, “Monaural speech segregation based on pitchtracking and amplitude modulation,” IEEE Trans. Neural Netw. , vol. 15,no. 5, pp. 1135–1150, 2004.[59] L. Atlas and C. Janssen, “Coherent modulation spectral filtering forsingle-channel music source separation,” in IEEE Proc. ICASSP’05 ,vol. 4, 2005, pp. 461–464. [60] P. X. Joris, C. E. Schreiner, and A. Rees, “Neural processing ofamplitude-modulated sounds,” Physiol. Rev. , vol. 84, no. 2, pp. 541–577, 2004.[61] F.-G. Zeng, K. Nie, G. S. Stickney, Y.-Y. Kong, M. Vongphoe, A. Bhar-gave, C. Wei, and K. Cao, “Speech recognition with amplitude andfrequency modulations,” PNAS , vol. 102, no. 7, pp. 2293–2298, 2005.[62] R. V. Shannon, F.-G. Zeng, V. Kamath, J. Wygonski, and M. Ekelid,“Speech recognition with primarily temporal cues,” Science , vol. 270,no. 5234, pp. 303–304, 1995.[63] B. R. Glasberg and B. C. J. Moore, “Derivation of auditory filter shapesfrom notched-noise data,” Hear. Res. , vol. 47, no. 1, pp. 103–138, 1990.[64] J. L. Flanagan, “Parametric coding of speech spectra,” J. Acoust. Soc.Am. , vol. 68, no. 2, pp. 412–419, 1980.[65] J. E. Shoup and L. L. Pfeifer, “Acoustic characteristics of speechsounds,” in Contemporary Issues in Experimental Phonetics . AcademicPress, 1976, pp. 171–224.[66] A. Keitel, J. Gross, and C. Kayser, “Perceptually relevant speech trackingin auditory and motor cortex reflects distinct linguistic features,” PLOSBiol. , vol. 16, no. 3, p. e2004473, 2018.[67] H. R. Bosker and M. Cooke, “Talkers produce more pronouncedamplitude modulations when speaking in noise,” J. Acoust. Soc. Am. ,vol. 143, no. 2, pp. EL121–EL126, 2018.[68] P. R. Hoskins, K. Martin, and A. Thrush, Diagnostic ultrasound: Physicsand equipment , 3rd ed. CRC Press, 2019.[69] M. Felsberg and G. Sommer, “The monogenic signal,” IEEE Trans.Signal Process. , vol. 49, no. 12, pp. 3136–3144, 2001.[70] H. Aragonda and C. S. Seelamantula, “Demodulation of narrowbandspeech spectrograms using the Riesz transform,” IEEE Audio, Speech,Language Process. , vol. 23, no. 11, pp. 1824–1834, 2015.[71] C. S. Seelamantula, N. Pavillon, C. Depeursinge, and M. Unser, “Localdemodulation of holograms using the Riesz transform with applicationto microscopy,” J. Opt. Soc. Am. A , vol. 29, no. 10, pp. 2118–2129,2012.[72] K. Nadeau, A. J. Durkin, and B. J. Tromberg, “Advanced demodulationtechnique for the extraction of tissue optical properties and structuralorientation contrast in the spatial frequency domain,” JBO , vol. 19, no. 5,p. 056013, 2014.[73] L. M. Sanchez-Brea and F. J. Torcal-Milla, “Near-field diffraction ofgratings with surface defects,” Appl. Opt. , vol. 49, no. 11, pp. 2190–2197, 2010. Supplementary Material for “Amplitude Demodulation of Wideband Signals” Mantas Gabrielaitis mantas.gabrielaitis@ { ist.ac.at, gmail.com } C ONTENTS Overview RECOVERY CONDITIONS A Auxiliary Proofs B Recovery Proofs C Further Analysis: Numerical Experiments AP ALGORITHMS D Mathematical Preliminaries E Properties of the Constraint Sets and Associated Metric Projections F Convergence Proofs G Lower Bound on the Convergence Error SIGNAL CLASSES H Types of Wideband Carriers Found in Practice I Synthetic Test Signals a r X i v : . [ ee ss . SP ] F e b PERFORMANCE TESTS J Configurations of Demodulation Algorithms for Performance Analysis K Implementation on a Computer ADDITIONAL RESULTS L Errors of Carrier Estimates M Convergence Rates at Different Signal Lengths N Repeated Demodulation O Demodulation of Speech Signals References VERVIEW 3 O VERVIEW This document serves as a source of additional information to support the ideas and results introduced in theaccompanying paper “Amplitude Demodulation of Wideband Signals.” Virtually, the material provided in this supplement can be divided into five blocks comprising, respectively, Sections A – C , D – G , H – I , J – K , and L – O : ‚ The A – C block provides proofs of modulator recovery conditions. ‚ The D – G block is concerned with mathematical aspects of the alternating projection (AP) algorithms ofamplitude demodulation. ‚ The H – I block reviews the main types of amplitude-modulated wideband signals found in practice anddefines the synthetic modulators and carriers used for testing purposes in the present work. ‚ The J – K block presents details on the numerical implementation of the AP and other relevant demodulationmethods on a computer, as well as their benchmarking configurations. ‚ The L – O block contains auxiliary simulation results and their discussion.A summary of each of the sections follows next to ease navigation through this document. Section A establishes several auxiliary results that are exploited in the modulator recovery proofs next. Section B provides proofs for the modulator recovery conditions introduced in Section II-C of the main text. Section C discusses the results and implementation of the numerical experiments performed to extend themodulator recovery conditions. Section D introduces the basic concepts of mathematical analysis necessary for the formulation and study ofthe properties of AP algorithms. Section E formulates and proves relevant properties of the constraint sets of modulator estimates and definesoperators that implement metric projections onto the modulator constraint sets used in this work. Section F provides the convergence proofs of the AP algorithms introduced in the main paper. Section G derives a lower bound on the convergence error. Section H reviews the main types of amplitude-modulated wideband signals found in practice. Section I defines the modulators and carriers of synthetic signals used to test amplitude demodulationalgorithms in the present work. Section J lists configurations of the execution control parameters employed for the performance analysis ofthe AP and LDC algorithms of amplitude demodulation. Section K provides information on the implementation and execution of the demodulation algorithms on acomputer. OVERVIEW Section L overviews the results of demodulation algorithm performance tests in terms of carrier estimates. Section M presents additional simulation results on the dependence of convergence rates of the AP algorithmson the signal length. Section N introduces simulation results of repetitive demodulation of carrier and modulator estimates obtainedusing the AP-A algorithm. Section O presents the results of additional simulations that demonstrate the suitability and consistency ofthe AP approach to demodulate wideband speech signals. ECOVERY CONDITIONS 5 RECOVERY CONDITIONS A. A UXILIARY P ROOFS Here, we establish two important properties of the unitary DFT basis vectors that are repetitively used in theproofs of propositions about the modulator recovery conditions in the next section. Proposition A.1. Consider a subset of DFT basis vectors t f p k q u k P I ω ˚ n . Assume a set of arbitrarily chosen n s sample points encoded by components of a vector r P N n s ` , and introduce a linear transform L r that mapsevery x P R n to r x r , x r , . . . , x r ns s T P R n s . Then, a set t L r f p k q u k P I ω ˚ n is linearly independent if and only if n s ě ω ˚ ´ .Proof. Necessity. If n s ă ω ˚ ´ , the set of vectors t L r f p k q u k P I ω ˚ n is linearly dependent because the number oflinearly independent vectors cannot be higher than the number of components they consists of.Sufficiency. First, consider the case of n s “ ω ˚ ´ . Then, t L r f p k q u k P I ω ˚ n is linearly independent if and onlyif the determinant of a matrix formed by concatenating all vectors from this set in an arbitrary order is not equalto zero (see, e.g., [1, p. 13]). To show that the latter condition is satisfied in our case, define a matrix M “ L r r f p n ´ ω ˚ ` q , f p n ´ ω ˚ ` q , . . . , f p n q , f p q , f p q , . . . , f p ω ˚ ´ q , f p ω ˚ q s . (24)Taking into account that L r f p k q can be written as rp z r q k ´ , p z r q k ´ , . . . , p z r ns q k ´ s T {? n , with z “ e ı π { n , andthat p z r q k “ p z r q k mod n , M can be expressed as a product of a diagonal matrix and a Vandermonde matrix: M “ diag r? n ¨ L r f p n ´ ω ˚ ` q s r L r f p q , L r f p q , . . . , L r f p ω ˚ ´ q s . (25)Thus, det M “ ? n ¨ det diag r L r f p n ´ ω ˚ ` q s ¨ det r L r f p q , L r f p q , . . . , L r f p ω ˚ ´ q s“ ? n ¨ ω ˚ ´ ź i “ ` L r f p n ´ ω ˚ ` q ˘ i ¨ ω ˚ ´ ź i “ i ´ ź j “ `` L r f p q q i ´ p L r f p q ˘ j ˘ “ ? n ¨ ω ˚ ´ ź i “ e ı π p n ´ ω ˚ ` q r i { n loooooooomoooooooon ‰ ¨ ω ˚ ´ ź i “ i ´ ź j “ ` e ı πr i { n ´ e ı πr j { n ˘looooooooooomooooooooooon ‰ ‰ , (26)which implies that the set t L r f p k q u k P I ω ˚ n is linearly independent for n s “ ω ˚ ´ . When writing the secondequality above, we used the expression of the determinant of a Vandermonde matrix (see, e.g., [1, p. 143]).It follows from the definition of linear independence (see, e.g., [1, p. 3]) that extending each vector in theset by additional components cannot change the set from linearly independent to linearly dependent. Hence, thelinear independence of t L r f p k q u k P I ω ˚ n for n s “ ω ˚ ´ implies the linear independence of t L r f p k q u k P I ω ˚ n for n s ą ω ˚ ´ . (cid:4) RECOVERY CONDITIONS Remark. The fact that t L r f p k q u k P I ω ˚ n is linearly independent for n s ě ω ˚ ´ means that the system of linearequations L r x “ ř k P I ω ˚ n ` a k ¨ L r f p k q ˘ has a unique solution for all x P S ω ˚ if n s ě ω ˚ ´ . In other words,every x P S ω ˚ can be recovered from its known sample points then. Proposition A.2. Consider some x p q P S ω ˚ and x p q P S ω ˚ . Assume a set of regularly chosen n s ě ω ˚ ´ sample points encoded by components of a vector r P N n s ` such that r i ` ´ r i “ n { n s for every i P I n s . Then, } x p q } {} x p q } “ } L r x p q } {} L r x p q } . (27) Proof. By the definition of S ω ˚ (see Section II-A in the main text), x “ ÿ k P I ω ˚ n ` a k ¨ f p k q ˘ , x P S ω ˚ . (28)Moreover, it follows from the unitary property of the DFT matrix, x f p j q , f p k q y “ δ j,k , that } x } “ ÿ k P I ω ˚ n | a k | . (29)Applying L r to both sides of (28) yields L r x “ ÿ k P I ω ˚ n ` a k ¨ L r f p k q ˘ , x P S ω ˚ , (30)where, taking into account that r i ` ´ r i “ n { n s , L r f p k q “ r e ı π ¨ r ¨p k ´ q{ n , e ı π ¨ r ¨p k ´ q{ n , . . . , e ı π ¨ r ns ¨p k ´ q{ n s T {? n “ r , e ı π ¨ ¨p k ´ q{ n s , . . . , e ı π ¨p n s ´ q¨p k ´ q{ n s s T {? n s ¨ ` e ı π ¨ r ¨p k ´ q{ n ¨ a n s { n ˘ . (31)Note that L r f p k q is the k -th column of the unitary n s ˆ n s DFT matrix multiplied by a coefficient whose absolutevalue is equal to a n s { n . Therefore, analogously to (29), } L r x } “ p n s { n q ¨ ÿ k P I ω ˚ n | a k | , (32)as long as n s ě ω ˚ ´ . Consequently, } L r x } “ a n s { n ¨ } x } for x P S ω ˚ , which implies (27). (cid:4) B. R ECOVERY P ROOFS In this section, we prove the modulator recovery conditions stated in the main text in the form of Proposi-tions II.1 – II.4 . We repeat the original assertions from the main text for the sake of convenience. Here, and in the sequel, δ i,j denotes the Kronecker delta. If n s ă ω ˚ ´ , some of the vectors L r f p k q and L r f p j q are identical for k ‰ j , and hence, (32) does not apply. ECOVERY CONDITIONS 7 Proposition II.1. For almost every m P M ω , ˆm “ m only if c P C d , $ ě ω , and n s ” ř ni “ I t u p| c i |q ě $ ` ω ´ . Proof. We prove the proposition by showing that, almost everywhere in M ω , m ‰ ˆm if c R C d , or $ ă ω , or n s ă $ ` ω ´ . For the sake of convenience, we restate the definition of ˆm here: ˆm “ arg min x P S ě| s | X S $ } x } . (33)If c R C d , it means that either | c i | ă , for every i P I n , or there exists at least one i P I n such that | c i | ą .In the former case, | s i |{ m i ă for every i P I n . Hence, for α “ max t| s i |{ m i u i P I n , α ¨ m belongs to S ě| s | X S $ but has a smaller norm than m , i.e., m ‰ ˆm . In the latter case, ˆ arg min x P S ě| s | X S $ } x } ˙ i ą m i (34)for all i corresponding to | c i | ą because | s i | ą m i then. Thus, m “ ˆm does not apply either. Therefore, m “ ˆm holds only if c P C d .If $ ă ω , then the subset of modulators for which m “ ˆm is valid has the cardinality of R $ ´ , and hence,has zero volume in M ω , whose cardinality is that of R ω ´ .Next, assume that c P C d , and $ ě ω , but n s ă $ ` ω ´ . Let us represent indexes of all sample pointscorresponding to | c i | “ by a vector r P N n s ` . Then, analogously to (30), we have L r m “ ÿ k P I $n ` a k ¨ L r f p k q ˘ . (35)For the sake of convenience, let us redefine f p k q “ $’’’’&’’’’% ϕ p q , k “ p ϕ p k ´ q ` ı ¨ ϕ p k ´ q q{? , ď k ď $ p ϕ p p n ´ k q` q ` ı ¨ ϕ p p n ´ k q` q q{? , n ´ $ ` ď k “ n , (36)and a k “ $’’’’&’’’’% α , k “ p α k ´ ` ı ¨ α k ´ q{? , ď k ď $ p α p n ´ k q` ` ı ¨ α p n ´ k q` q{? , n ´ $ ` ď k “ n . (37)Then, (35) turns into L r m “ $ ´ ÿ i “ ` α i ¨ L r ϕ p i q ˘ , (38)According to Proposition A.1 , a set of n s vectors t L r f p q , . . . , L r f p r p n s ` q{ s , L r f p n ´ t p n s ´ q{ u , . . . , L r f n u islinearly independent. Hence, by (36), the same applies to t L r ϕ p i q u n s i “ . Consequently, (38) is an underdetermined Here ˆm is as defined by (6) in the main text. RECOVERY CONDITIONS system of linear equations defined by a full-rank matrix. We know from linear algebra that a general solutionof such system is expressed as a sum of its separate solution α p q and a solution of “ $ ´ ÿ i “ ` α i ¨ L r ϕ i ˘ , (39)Solutions of (39) form a p $ ´ ´ n s q -dimensional subspace of R $ ´ . Thus, we can express the generalsolution of (38) as α “ α p q ` $ ´ ´ n s ÿ i “ z i ρ p i q , z P R $ ´ ´ n s , (40)where t ρ p i q u $ ´ ´ n s i “ is an orthonormal basis of the space of solutions of (39). Taking into account (36) – (37)as well as the linear independence of t f p k q u k P I $n and t ρ p i q u $ ´ ´ n s i “ , (40) together with (30) define a linearinjective function that maps from R $ ´ ´ n s to S $ : f p z q “ $ ´ ÿ j “ ˆ α p q ` $ ´ ´ n s ÿ i “ z i ρ p i q ˙ j ϕ p j q , z P R $ ´ ´ n s . (41)The image of f p z q is a subset of those elements of S $ that coincide with the true modulator m at entries r P N n s ` . The injective nature of this function guarantees the existence of a unique z m P R $ ´ ´ n s such that m “ f p z m q .Now, assume an m P M ` ω “ t x P M ω : m i ą , i P I n u , i.e., any feasible modulator whose all entriesare strictly positive. Define an (cid:15) “ min t m i ´ | s i | : p| c i | ă q ^ p i P I n qu . (cid:15) exists and is positive because | c i | “ only for n s ă n out of n components of c by the assumption of the proposition, and m i ´ | s i | “ m i p ´ c i q . The linearity of f p z q implies its continuity at every point of its domain. Hence, there exists an η such that } f p z q ´ f p z m q} ă (cid:15) whenever z P H z m ,η “ t z P R $ ´ ´ n s : } z ´ z m } ă η u . On the other hand, } f p z q ´ f p z m q} “ ař ni “ rp f p z qq i ´ m i s ă (cid:15) implies |p f p z qq i ´ m i | ă (cid:15) , and hence p f p z qq i ą | s i | , for every i P I n . Consequently, f p z q P p S ě| s | X S $ q , z P H z m ,η . (42)Furthermore, if m “ ˆm , then } m } ă } f p z q} for every z P p H z m ,η z z m q , i.e., z m is a strict local minimumpoint of } f p z q} “ ››››› $ ´ ÿ j “ ˆ α p q ` $ ´ ´ n s ÿ i “ z i ρ p i q ˙ j ϕ p j q ››››› “ ››››› $ ´ ÿ j “ ˆ α p q ` $ ´ ´ n s ÿ i “ z i ρ p i q ˙››››› “ } α p q } ` $ ´ ´ n s ÿ i “ ` z i ` ¨ z i ¨ x α p q , ρ p i q y ˘ . (43) ECOVERY CONDITIONS 9 } f p z q} is a continuous, differentiable function with a positive-definite Hessian: B } f p z q} {B z i B z j “ δ i,j . Thus,it has a unique strict local minimum point defined by B} f p z q} {B z i “ : z : i “ ´x α p q , ρ p i q y , i P I $ ´ ´ n s . (44)Remember that, without loss of generality, α p q is a particular solution of (38). α p q corresponding to m , i.e., α p q ” p α p q , α p q , . . . , α ω ´ , , . . . , q T with m “ ř ω ´ i “ α p q i ϕ p i q , is exactly such a solution. In this case, asfollows from (41), m “ f p z : q if and only if z : “ . According to (44), that is equivalent to requiring α p q tobe a solution of the following homogeneous system of linear equations: x α p q , ρ p i q y “ , i P I $ ´ ´ n s . (45)Hence, the subset of M ` ω to which m “ ˆm applies has the same cardinality as R D , where D is the dimensionof the solution space of (45). Taking into account the linear independence of t ρ p i q u $ ´ ´ n s i “ , D is equal tothe difference between the number of elements of α p q that are not identically equal to zero and the numberof equations that are not trivially satisfied by any feasible α p q . The latter depends on c , specifically, on thepositions of sample points with | c i | “ . To see this, consider two cases: ‚ A carrier with equidistantly-spaced true sample points: | c i `p j ´ q¨ d | “ for some i P I d and every j P I n s ,where d “ n { n s . Then, it follows from (31) that some of the elements of the system t L r f p k q u k P I $n areidentical as long as n s ă $ ´ . Specifically, we have L r f p k q “ $&% L r f p k ` n ´ n s q , t p n s ` q{ u ď k ď $ L r f p k ` n s ´ n q , n ´ $ ` ď k ď n ´ t p n s ´ q{ u . (46)Equivalently, L r ϕ n s ` χ ns ´ i “ p´ q χ i ¨ L r ϕ n s ´ χ ns ` χ i ` i , ď i ď $ ´ ` χ n s ´ n s (47)and L r ϕ p n s ` q “ when χ n s “ , where, χ z “ z mod . Consequently, p ρ p i q q j “ $’’’’&’’’’% {? , j “ n s ` χ n s ´ i p´ q χ i ` {? , j “ n s ´ χ n s ` χ i ` i , otherwise , ď i ď $ ´ ` χ n s ´ n s , (48)and ρ p $ ´ ´ n s q “ when χ n s “ . Moreover, by our choice of α p q , α p q j “ , ω ď j ď $ ´ . (49) This strict local minimum point is the only local minimum point of this function. Note that t ρ i u $ ´ ´ n s i “ is linearly independent. Note that n s ă $ ` ω ´ and $ ě ω imply n s ă $ ´ . Note that, as discussed before, t ϕ i u n s i “ is linearly independent. (48) and (49) together imply that (45) holds for i P I n s ´p ω ´ q independent of α p q corresponding to thechosen m P M ` ω , and that (45) holds for the remaining i P p I $ ´ ´ n s z I n s ´p ω ´ q q if and only if α p q j “ , p n s ` ´ $ q ď j ď ω ´ . (50)Hence, (45) applies only to the subset of M ` ω with D “ min tp ω ´ q , p ω ´ q ´ p $ ` ω ´ ´ n s qu“ min tp ω ´ q , p n s ` ´ $ q ´ u . (51)(51) implies that D ă ω ´ as long as n s ă $ ` ω ´ . Then, the subset of all m P M ` ω that satisfy m “ ˆm has zero volume in M ` ω , which has the cardinality of R ω ´ . Moreover, if n s ă $ , it followsfrom (50) that (45) has only the trivial solution α p q “ , which implies m “ . The latter is infeasiblein our case, and hence, m “ ˆm does not apply to any m P M ` ω . On the other hand, if n s ě $ ` ω ´ ,then all equations of the (45) system are satisfied independent of the actual feasible α p q . Consequently, D “ ω ´ . ‚ An arbitrary c P C d with unspecified structure. Then, none of the $ ´ ´ n s equations of the system (45)are satisfied independent of the actual α p q . Therefore, D “ p ω ´ q ´ p $ ´ ´ n s q“ n s ´ p $ ´ ω q . (52)(52) is valid only if p ω ´ q ą p $ ´ ´ n s q . Otherwise, (45) has only the trivial solution α p q “ , whichimplies m “ . The latter is infeasible in our case, and hence, m “ ˆm does not apply to any m P M ` ω . Inany case, D ă p ω ´ q as long as n s ă p $ ´ q , which follows from the assumption of the propositionthat n s ă p $ ` ω ´ q and $ ě ω . In fact, p $ ´ q ´ p $ ` ω ´ q “ p $ ´ ω q .Compared with the general case, a smaller number of necessary sample points with | c i | “ is achieved inthe first example due to the fact that the subset of vectors L r ϕ p k q with ω ď k ď $ ´ is linearly dependentand that α p q k are identically equal to zero in that range. Specifically, every ϕ p k q in the range ω ď n s ´ ` χ n s is proportional to one of the ϕ p k q in the range n s ` ´ χ n s ď $ ´ . Every such dependence reduces thenecessary number of points with | c i | “ for modulator recovery by one. Further, it follows from (31), (36),and Proposition A.1 that the sets t L r ϕ p k q u n s ´ ` χ ns k “ ω and t L r ϕ p k q u $ ´ k “ n s ` ´ χ ns are linearly independent. Hence,no further linear dependencies among t L r ϕ p k q u $ ´ k “ ω are possible in general. This means that n s “ $ ` ω ´ is the absolute minimum of sample points with | c i | “ necessary for m “ ˆm to hold. The second exampleabove illustrates that this number is surely higher for some c P C d .The last three paragraphs demonstrate that the subset of M ` ω to which m “ ˆm applies has zero volume in M ` ω if n s ă $ ` ω ´ . Taking into account that p M ω z M ` ω q has lower cardinality than M ` ω ( R ω ´ vs. R ω ´ ),we conclude that the set of all m P M ω that satisfy m “ ˆm has zero volume in M ω if n s ă $ ` ω ´ . One possible example of this kind is c with | c i | “ for i P I n s , and | c i | ă for i P p I n z I n s q . ECOVERY CONDITIONS 11 (cid:4) Proposition II.2. Consider m P M ω and ˜c P C d with | ˜ c i | “ for i P J n Ď I n . If ˆm “ m holds for the m and ˜c , then it also holds for every pair made of the same m and any c P C d with | c i | “ for i P J n .Proof. Denote ˜s “ m ˝ ˜c and s “ m ˝ c . Now, note that | c i | ě | ˜ c i | by the condition of the proposition,and hence, | s i | ě | ˜ s i | . Therefore, S ě| s | Ď S ě| ˜s | , which implies p S ě| s | X S $ q Ď p S ě| ˜s | X S $ q , and, consequently, arg min x P S ě| s | X S $ } x } ě arg min x P S ě| ˜s | X S $ } x } . According to the proposition, arg min x P S ě| ˜s | X S $ } x } “ m . On the other hand, byconstruction, m P p S ě| s | X S $ q . Thus, arg min x P S ě| s | X S $ } x } “ m . (cid:4) Remark. It can be shown by example that the validity of ˆm “ m p q for some m p q P M ω and ˜c P C d with | ˜ c i | “ , i P J n , does not necessarily imply the validity of ˆm “ m p q for another m p q P M ω and the same ˜c . Proposition II.3. Assume m P M ω and c P C d with $ ě ω . If, additionally, there exist d P I n and i P I d suchthat n s ” p n { d q P N ` , n s ě $ ` ω ´ , and | c i `p j ´ q¨ d | “ for every j P I n s , then ˆm “ m . Proof. Here, we distinguish between two cases: p $ ` ω ´ q ď n s ă p $ ´ q and n s ě p $ ´ q .If p $ ` ω ´ q ď n s ă p $ ´ q , then it follows from the proof of Proposition II.1 that m has the smallestnorm among all elements of the image of f p z q defined by (41), i.e, all elements x P S $ that satisfy L r x “ L r m .Next, consider a y P S ω such that p L r y q i ě p L r m q i for every i P I n s and p L r y q i ą p L r m q i for at least one i P I n s . Then, by (27), } y } ą } m } . Moreover, using the same argumentation as for m , we see that y hasthe smallest norm among all elements x P S $ that satisfy L r x “ L r y . Hence, m has smaller norm than anyother element x P S $ that satisfies L r x ě L r m . Moreover, p S ě| s | X S $ q Ă S $ . Thus, we conclude that ˆm “ m holds.If n s ě p $ ´ q , then it follows from (27) of Proposition A.2 that arg min x P S ě| s | X S $ } x } “ arg min x P S ě| s | X S $ } L r x } . (53)Further, the constraint set S ě| s | implies that } L r x } ě } L r | s |} , x P S ě| s | X S $ . (54)On the other hand, | s r i | “ | m r i ¨ c r i | “ m r i ¨ | c r i | “ m r i for i P I n s , i.e., L r | s | “ r m r , m r , . . . , m r ns s T .According to Proposition A.1 , t L r f p k q u k P I $n is linearly independent if n s ě $ ´ . Therefore, (30) has aunique solution, which, by (28), means that, among all S $ , there is a unique x “ m that satisfies L r x “r m r , m r , . . . , m r ns s T “ L r | s | . Hence, by (54), arg min x P S ě| s | X S $ } L r x } “ m , and, by (53), arg min x P S ě| s | X S $ } x } “ m . (cid:4) Here ˆm is as defined by (6) in the main text. Note that ω ˚ in (28) – (30) stands for $ here. Proposition II.4. Consider m P M ω and c P S | .. |ď . Take n s ě $ ´ sample points of s “ m ˝ c whoseindexes are defined as entries of any chosen r P N n s ` with r i ` ´ r i “ n { n s for every i P I n s . Then, } m ´ ˆm } {} m } ď b ´ ř n s i “ s r i { ř n s i “ m r i . (55) Proof. For the sake of convenience, we will exploit the linear transformation L r , already introduced in Propo-sition A.1 , that maps every x P R n to r x r , x r , . . . , x r ns s T . Then, (55) can be rewritten as } m ´ ˆm } {} m } ď b ´ } L r s } {} L r m } . (56)Note that } m } ´ } ˆm } ´ } m ´ ˆm } “ ¨ } ˆm } ¨ p} m } ´ } ˆm } q . (57)Next, we have by construction that m P p S ě| s | X S $ q . Hence, } m } ě } ˆm } , (58)which, together with (57) implies } m } ´ } ˆm } ´ } m ´ ˆm } ě , i.e., } m ´ ˆm } {} m } ď b ´ } ˆm } {} m } . (59)On the other hand, by Proposition A.2 , } ˆm } {} m } “ } L r ˆm } {} L r m } if n s ě $ ´ and r i ` ´ r i “ n { n s for every i P I n s . Thus, } m ´ ˆm } {} m } ď b ´ } L r ˆm } {} L r m } . (60)Finally, ˆ m i ě | s i | for every i P I n because m P S ě| s | , which means that } L r ˆm } ě } L r s } . (61)Combining (60) with (61) leads to (56). (cid:4) Remark. Note that } m } “ } ˆm } in (58) if and only if m “ ˆm because S ě| s | X S $ is convex and } . . . } isstrictly convex. Therefore, the equality in (55) holds if and only if m “ ˆm , i.e., the modulator recovery is exact.C. F URTHER A NALYSIS : N UMERICAL E XPERIMENTS Here, we present the results of numerical experiments used to extend the modulator recovery conditions tocarriers with nonuniformly placed sample points | c i | “ in terms of the parameters n , ω , $ , and d . Here ˆm is as defined by (6) in the main text. ECOVERY CONDITIONS 13 Setup The numerical experiments under consideration consist of the following steps.1. pairs of m and c are generated by randomly sampling from M ω and C d for every feasible combinationof ω and d consistent with a chosen signal length n .2. For every pair of m and c generated, m is inferred from the s “ m ˝ c via ˆm defined by (6). The latter isevaluated by using the AP-P algorithm, introduced in Section III-C of the main paper, with (cid:15) tol “ ´ and unlimited N iter .3. For every combination of the parameters ω and d , two estimates related to the recovery error are evaluated:1) the average empirical error x E m y and 2) the fraction of cases with vanishing error P p E m ă ε q , where ε is a positive number arbitrarily close to zero. P p E m ă ε q can be seen as the demodulation success ratefor a given error threshold ε .A crucial aspect of the above experiments to producing informative data for our purposes is the way M ω and C d are sampled. For both of these sets, we exploited uniform sampling but with some additional constraints, asexplained next. ‚ The cutoff frequency ω , defining the modulator set, and the cutoff frequency $ , defining the estimator ˆm ,were fixed to be equal. This choice allowed us to considerably reduce the extent of relevant parametercombinations to be checked without loss of generality. Indeed, $ ě ω is a necessary condition for a fullrecovery independent of c P C d and M ω Ă M $ if $ ą ω . Hence, all recovery conditions applicable inthe case of $ “ ω hold for $ ą ω as well. ‚ Only the subset of pure spike-train carriers consisting of c i P t , u sample points was considered amongall possible c P C d . According to Proposition II.2 , that is sufficient for identifying full recovery conditionswithout loss of generality. ‚ Different elements of the pure spike-train subset of C d may substantially differ in the number n s of samplepoints with | c i | “ . In particular, r n { d s ď n s ď n ´ d ` . We considered modulator reconstructionby uniformly sampling either from the sparsest ( n s “ r n { d s ) or the densest ( n s “ n ´ d ` ) subset ofspike-train carriers. In view of Proposition II.2 , pure spike-train carriers with n s “ r n { d s have the tightest,and hence the most general, constraints for exact modulator recovery in terms of the parameters $ and d .The algorithms for sampling from the modulator and carrier sets specified above are presented in Section I. Results Fig. 10 A displays color-plots of the fraction P of the modulator recovery cases with E m ă ´ overthe p d, $ q plane for n “ and n s “ r n { d s . In agreement with the necessary recovery condition discussedin Section II-C, P p E m ă ´ q is equal to 0 for all points with r n { d s ă $ ´ . More remarkably, for r n { d s ě $ ´ , P p E m ă ´ q equals 1 except some boundary points, where P varies between 0.65 and 1. n = 32 n = 256 d d d d d d d d 〈 E m 〉 n s = ⌈ n/d ⌉ P(E m < -12 ) n s = ⌈ n/d ⌉ P(E m < -12 ) n s = ⌈ n/d ⌉ 〈 E m 〉 n s = ⌈ n/d ⌉ C B A D 〈 E m 〉 n s =n − d+1 P(E m < -12 ) n s =n − d+1 P(E m < -12 ) n s =n − d+1 〈 E m 〉 n s =n − d+1 C B A D Fig. 10. Success rates and errors of modulator recovery for signals based on pure c i P t , u spike-train carriers. A – A : color plots ofthe fraction ( P ) of modulator recovery cases with the recovery error E m lower than ´ for different combinations of d and $ , and n “ ; red lines plot the relation r n { d s “ $ ´ . A displays the results for carriers with n s “ r n { d s spikes, while A correspondsto carriers with n s “ n ´ d ` . B – B : the same as A – A , but with the average recovery error instead of the success rate over allmodulator and carrier pairs shown for each combination of d and $ . C – C : the same as A – A except that n “ . D – D : thesame as B – B except that n “ . However, even in the latter cases, the error is small, as follows from Fig. 10 B , which plots the average x E m y over the p d, ω q plane. The maximum likelihood (ML) estimate of P p E m ă ´ q for all tested modulator-carrierpairs that adhere to the necessary recovery condition r n { d s ě $ ´ is 0.971; its confidence interval is p . , . q .Increasing the number of carrier points with | c i | “ does not change the landscape of the recovery successrate considerably. Indeed, pushing n s to the maximum n ´ d ` increases the P values only at points in theimmediate vicinity of the r n { d s “ $ ´ boundary (see Fig. 10 A ). Nevertheless, it has to be noted that therecovery errors are decreased by increasing n s on average (see Fig. 10 B ). The ML estimate of P p E m ă ´ q for all tested modulator-carrier pairs that adhere to the necessary recovery condition r n { d s ě $ ´ is 0.987in this case; its confidence interval is p . , . q .We found an analogous picture while considering signals with different lengths n . One such example ( n “ )is considered in Fig. 10 C – D . The chances of the modulator recovery with vanishing error, i.e., E m ă ´ ,are higher in this case. Specifically, the ML estimate of P p E m ă ´ q is 0.9933, with the confidenceinterval p . , . q when n s “ r n { d s . That can be explained by the smaller contribution of the boundary ECOVERY CONDITIONS 15 points of the relation r n { d s “ $ ´ in the p d, $ q plane to the total count. It is important to note that, in allcases discussed here, essentially the same results are obtained even if the error threshold ε is increased to ´ .This rejects any possibility of numerical inaccuracies affecting our conclusion. AP ALGORITHMS D. M ATHEMATICAL P RELIMINARIES In this section, we introduce some basic concepts of mathematical analysis necessary for the formulation andassessment of the AP algorithms that we used to calculate modulator estimators defined by (6) and (8) in themain paper. Convex, interior, closed, and bounded sets We start with definitions of a few basic attributes of sets in Euclidean spaces. Definition D.1. A set S Ď R n is said to be convex if θ ¨ x ` p ´ θ q ¨ y P S for all x , y P S and θ P r , s . Definition D.2. An element x P S Ď R n is said to be an interior point of S if there exists an (cid:15) ą such that t y P R n : } x ´ y } ă (cid:15) u Ă S . Definition D.3. A set that consists of all interior points of S Ď R n is called the interior of S . We denote it by S ˝ . Definition D.4. An element y P R n is said to be a contact point of S Ď R n if, for any (cid:15) ą , there exists an x P S such that } x ´ y } ă (cid:15) . Definition D.5. A set S Ď R n is said to be closed if it is equal to the set of all its contact points.Convexity and closedness of sets are preserved under intersection. Proposition D.6. The intersection S X S of two closed and convex sets S and S is closed and convex. Another important characteristic of sets is their boundedness. Definition D.7. A set S Ă R n is said to be bounded if there exists a b P R such that } x ´ y } ď b for all x , y P S . Definition D.8. A set S Ă R “ R is said to be bounded from above if there exists a u P R such that x ď u for all x P S . u is called an upper bound of S . Remark. u P R is said to be the least upper bound of S Ă R if x ď u for all x P S , and there exists a y P S for every (cid:15) ą such that y ą u ´ (cid:15) . The least upper bound exists for any S Ă R bounded from above due tothe continuity of the real numbers. Definition D.9. A set S Ă R “ R is said to be bounded from below if there exists an l P R such that l ď x for all x P S . l is called a lower bound of S . P ALGORITHMS 17 Remark. ¯ l P R is said to be the greatest lower bound of S Ă R if ¯ l ď x for all x P S , and there exists a y P S for every (cid:15) ą such that y ă ¯ l ` (cid:15) . Analogously to the least upper bound, any S Ă R bounded from belowhas the greatest lower bound. Convergence of sequences The following fundamental properties of infinite sequences of points in bounded subsets of Euclidean spacesplay a critical role in the proofs of the convergence of the AP algorithms. Proposition D.10 (Monotone Convergence Theorem) . Any monotonically decreasing sequence of real numbers x p q ě x p q ě . . . ě x p i q ě . . . (62) that is bounded from below converges to its greatest lower bound ¯ l , i.e., for every (cid:15) ą , there exists an N p (cid:15) q such that | x p i q ´ ¯ l | ă (cid:15) whenever i ą N p (cid:15) q .Remark. Analogously, any monotonically increasing sequence that is bounded from above converges to its leastupper bound. Definition D.11. Consider a sequence x p q , x p q , . . . , x p i q , . . . in R n , i.e., x p i q P R n for every i ě . Anothersequence x p k q , x p k q , . . . , x p k i q , . . . in R n generated by removing some of the elements of the original sequenceis called a subsequence of the latter. Note that k i ą k j for all i ą j ě , and k i ě i for every i ě here. Proposition D.12 (Bolzano-Weierstrass Theorem) . Any bounded infinite sequence x p q , x p q , . . . , x p i q , . . . in R n has an infinite subsequence x p k q , x p k q , . . . , x p k i q , . . . that converges to a particular x : P R n , i.e., for any (cid:15) ą ,there exists an N p (cid:15) q such that } x p k i q ´ x : } ă (cid:15) whenever i ą N p (cid:15) q .Metric projections The central operation around which AP algorithms are built is that of a metric projection. Definition D.13. An element x z of a closed subset S of R n is said to be a metric projection of z P R n onto S if } x z ´ z } ď } x ´ z } for all x P S . We denote a transformation that assigns an x z P S to every z P R n by P S : R n Ñ S . Remark. P S generalizes a linear projection operator that assigns an element of a linear space to one of itssubspaces (see, e.g., [2, p. 223]). For the sake of brevity, we skip the qualifier “metric” and refer to P S as “aprojection” in the sequel. Proposition D.14 (see Theoreom 5.11 in [3]) . A projection of any element of R n onto its closed convex subset S exists and is unique. Inequalities Two important inequalities that we use extensively in the convergence proofs of AP algorithms apply toEuclidean spaces. Proposition D.15 (Triangle Inequality) . } x ` y } ď } x } ` } y } @ x P R n , @ y P R n . (63) Remark. Another inequality relevant to us follows from (63). In particular, let us consider some a , b , c P R n . If wedefine x “ a ´ b and y “ b ´ c , then (63) implies } a ´ c } ď } a ´ b } `} b ´ c } , i.e., } a ´ b } ě } c ´ a } ´} c ´ b } .On the other hand, setting x “ c ´ a and y “ a ´ b , we obtain } c ´ b } ď } c ´ a } ` } a ´ b } , i.e., } a ´ b } ě ´p} c ´ a } ´ } c ´ b } q . Thus, } a ´ b } ě ˇˇ } c ´ a } ´ } c ´ b } ˇˇ @ a P R n , @ b P R n , @ c P R n . (64) Proposition D.16 (Containing-Half-Space Inequality, see Theoreom 5.13 in [3]) . If S Ă R n is closed andconvex, then x x ´ P S r x s , P S r x s ´ y y ě @ x P R n , @ y P S . (65) Remark. S belongs to a half-space H x “ t z P R n : x x ´ P S r x s , P S r x s ´ z y ě u , which is defined for every x P p R n z S q . P S r x s is a boundary point of the H x .E. P ROPERTIES OF THE C ONSTRAINT S ETS AND A SSOCIATED M ETRIC P ROJECTIONS In this section, we establish convexity, closedness, and other relevant properties of the constraint sets S ě| s | and S $ that lay the basis for the formulation of the AP demodulation algorithms and determine their convergenceproperties. We then define the concrete metric projection operators of points in R n onto S ě| s | and S $ , which arethe main building blocks of the AP demodulation algorithms defined in Section III of the main paper. Properties of the constraint sets The constraint sets S ě| s | and S $ that define the AP approach to demodulation introduced in Section II.A ofthe main paper have the following properties. Proposition E.1. The set S ě| s | is convex and closed. Its interior is S ˝ ě| s | “ t x P R n : x i ą | s i | , i P I n u .Proof. The range of values of each component x i of x P S ě| s | and y i of y P S ě| s | is r| s i | , `8q Ă R . Obviously, θ ¨ x i ` p ´ θ q ¨ y i ě min r x i , y i s ě | s i | and θ ¨ x i ` p ´ θ q ¨ y i ď max r x i , y i s ă `8 , P ALGORITHMS 19 where θ P r , s . Therefore, θ ¨ x i ` p ´ θ q ¨ y i P r| s i | , `8q @ θ P r , s , @ i P I n . This, in turn, implies θ ¨ x ` p ´ θ q ¨ y P S ě| s | because the range of values of each component of x and y isindependent of the actual values of the remaining components. Hence, by Definition D.1 , S ě| s | is convex.To show that S ě| s | is closed, we first note that its complement S ě| s | “ R n z S ě| s | is equal to t x P R n : x i ă| s i | , i P I n u . For any x P S ě| s | and (cid:15) ď min r| s | ´ x s , there exists no y P S ě| s | such that } x ´ y } ă (cid:15) . Thus,none of the elements of S ě| s | are contact points of S ě| s | . On the other hand, trivially, all points of S ě| s | are itscontact points. Hence, S ě| s | coincides with the set of its contact points, i.e., it is a closed set.To determine the interior of S ě| s | , let us denote A “ t x P R n : x i ą | s i | , i P I n u , A i “ t x P R n : x i “ | s i | , x j ě | s j | , j P p I n zt i uqu @ i P I n . By construction, S ě| s | “ A Y pY ni “ A i q . Next, we note that t y P R n : } x ´ y } ă (cid:15) u Ă S ě| s | @ (cid:15) ď min r x ´ | s |s , @ x P A . Thus, by Definition D.2 , all points of A are interior points of S ě| s | . On the other hand, for all x P A i and (cid:15) ą , there exists a y P R n such that } x ´ y } ă (cid:15) and y R S ě| s | . In particular, this is satisfied by y suchthat y i “ x i ´ ¯ (cid:15) with ă ¯ (cid:15) ă (cid:15) and y j “ x j for all j P p I n zt i uq . Therefore, none of x P A i with i P I n areinterior points of S ě| s | . Altogether, this allows us to conclude that the interior of S ě| s | is equal to A . That A is nonempty follows straightforwardly from the fact that there exists an x i ą s i for any s i P R . (cid:4) Proposition E.2. The set S $ is convex, closed, and void of interior points. In particular, S $ is a linear subspaceof R n .Proof. It follows directly from the definition of R n that θ ¨ x ` p ´ θ q ¨ y P R n @ x P R n , @ y P R n , @ θ P r , s , (66)and thus, θ ¨ x $ ` p ´ θ q ¨ y $ P R n @ x ω P S $ , @ y $ P S $ , @ θ P r , s . (67)The definition of S $ implies that p Fx $ q i “ and p Fy $ q i “ for all x $ P S $ , y $ P S $ , and i P p I n z I $n q ,so that θ ¨ p Fx $ q i ` p ´ θ q ¨ p Fy $ q i “ p F p θ ¨ x $ ` p ´ θ q ¨ y $ qq i “ @ i P p I n z I $n q (68) Note that min r| s | ´ x s ą . Note that S $ Ă R n . because of the linearity of the Fourier transform. Combining (67) and (68) with the definition of S $ , we concludethat θ ¨ x $ ` p ´ θ q ¨ y $ P S $ for all x $ P S $ , y $ P S $ , i.e., S $ is convex.To prove that S $ is closed, we show that no point of S $ is a contact point of S $ . Indeed, the complementof S $ in R n is given by S $ “ y P R n : ř i Pp I n z I $n q p Fy q i ą ( . (69)Let us consider some y P S $ . Taking into account the definitions of S $ and S $ and the fact that F is unitary,we have that, for any x P S $ , } x ´ y } “ } F p x ´ y q} “ } Fx ´ Fy } “ ř i P I $n pp Fx q i ´ p Fy q i q ` ř i Pp I n z I $n q p Fy q i ě ř i Pp I n z I $n q p Fy q i ą . (70)Thus, for every y P S $ , there exist no x P S $ such that } x ´ y } ă (cid:15) with (cid:15) “ bř i Pp I n z I $n q p Fy q i , whichmeans that none of the elements of S $ are contact points of S $ . Therefore, S $ coincides with the set of itscontact points, i.e., it is closed.It follows from the definition of S $ that y “ p x ` (cid:15) ¨ e p q q R S $ for all x P S $ and (cid:15) ą , where e p q is theunit vector with all but its first components equal to zero. Indeed, p Fe p q q i “ {? n ‰ for all i P I n . Moreover,in that case, } x ´ y } “ (cid:15) . Thus, for all x P S $ and (cid:15) ą , there exists a y P R n such that } x ´ y } ă (cid:15) and y R S $ , which means that S ě| s | has no interior points.A necessary and sufficient condition for a subset S of a linear space R n to be a subspace is that p α ¨ x ` β ¨ y q P S for all α P R , β P R , x P S , and y P S (see, e.g., [2, p. 121]). That this applies to S $ follows from the proofof its convexity above if we replace θ and p ´ θ q by, respectively, α and β . (cid:4) Proposition E.3. The intersection of sets S ˝ ě| s | and S $ is nonempty, i.e., there exists an x P S ˝ ě| s | X S $ .Consequently, S ě| s | X S $ is also nonempty.Proof. Let us consider x “ λ ¨ , where is an element of R n with all its components equal to 1, and λ ą max r s s . It follows directly from the definition of S ě| s | that x P S ˝ ě| s | Ă S ě| s | . It also follows from thedefinitions of S $ and unitary discrete Fourier transform that p Fx q i “ δ i, ¨ ? n ¨ λ , i.e., x P S $ . Therefore, x P p S ˝ ě| s | X S $ q Ă p S ě| s | X S $ q . (cid:4) Remark. It is a direct consequence of Propositions D.6 , E.1 , and E.2 that S ě| s | X S $ is also closed and convex. Metric projections onto S ě| s | and S $ The metric projections of any point in R n onto its convex subsets relevant to us, i.e., S ě| s | and S $ , are achievedby the following operators. P ALGORITHMS 21 Proposition E.4. The metric projection operator P S ě| s | is defined by the elementwise maximum of the targetsignal s and the input argument z : P S ě| s | r z s “ | s | ` p z ´ | s |q ˝ θ p z ´ | s |q . (71) Proof. Note that ` P S ě| s | r z s ˘ i “ $&% z i , if z i ě | s i || s i | , if z i ă | s i | @ i P I n . (72)Hence, a necessary and sufficient condition for transforming any z P R n to x P S ě| s | is to increase everycomponent z i of z that does not satisfy z i ě | s i | by at least | s i | ´ z i , independent of values of the remainingcomponents. Now, we have from the definition of the Euclidean norm that } z ´ x } “ ař ni “ p z i ´ x i q @ z P R n , @ x P S ě| s | . (73)Thus, for every z P R n , } z ´ x } is minimized by an x P S ě| s | that is obtained by incrementing all components of z that satisfy z i ă | s i | by no more than necessary, i.e., by | s i | ´ z i , and leaving the remaining components intact.However, this is precisely how the operator P S ě| s | is defined via (72). Therefore, by using the Definition D.13 of the projection, we conclude that P S ě| s | projects z P R n onto S ě| s | . (cid:4) Proposition E.5. The metric projection operator P S $ is defined by a rectangular low-pass-filter transformation P S $ r z s “ p F ´ W $ F q z “ ÿ i P I $n x f p i q , z y ¨ f p i q . (74) Here, f p i q is the i-th column of the DFT matrix F . W $ is a diagonal matrix such that p W $ q ii “ $&% , if i P I $n , otherwise . (75) Proof. Let us consider an element of the set S $ expressed by x “ ř i P I $n a i ¨ f p i q . It follows from the definitionof the Euclidean norm (see Section II.E in the main paper) that } x ´ z } “ Bˆ z ´ ÿ i P I $n a i ¨ f p i q ˙ , ˆ z ´ ÿ i P I $n a i ¨ f p i q ˙F “ x z , z y ´ ¨ ÿ i P I $n a i ¨ x f p i q , z y ` ÿ i P I $n a i “ x z , z y ´ ÿ i P I $n x f p i q , z y ` ÿ i P I $n p a i ´ x f p i q , z yq . (76)When writing the second equality above, we used the fact that t f p q , f p q , . . . , f p n q u are orthonormal. It followsfrom the last equality of (76) that } x ´ z } , as a function of a i , is minimized by a i “ x f p i q , z y for every i P I $n .Thus, by Definition D.13 , ř i P I $n x f p i q , z y ¨ f p i q “ P S $ r z s is the projection of z onto S $ . (cid:4) Remark. P S $ is a linear operator, whereas P S ě| s | is not. F. C ONVERGENCE P ROOFS Here, we provide proofs of the propositions concerning the convergence of the AP algorithms that areformulated in Section III of the main paper. The proofs are adapted for finite-dimensional Euclidean spacesand exploit the particular structure of the modulator constraint sets. For the sake of convenience, we repeatthe original assertions as well. AP-B algorithm Proposition III.1. A sequence m p q , m p q , . . . , m p i q , . . . formed by the AP-B algorithm for (cid:15) tol “ and N iter Ñ`8 converges to some m : P S ě| s | X S $ . The convergence is geometric and monotonic, i.e., there exist γ ą and ă r ă such that } m p i q ´ m : } ď γ ¨ r i and } m p i ` q ´ m : } ď } m p i q ´ m : } for i ě .Proof. If the sequence m p q , m p q , . . . , m p i q , . . . terminates with some m p N q , i.e., m p N ` j q “ m p N q for every j ą , it follows from the formulation of the AP-B algorithm that m p N q “ a p N q , i.e., m p N q P S ě| s | X S $ . Thus, m p N q “ m : , which means that the solution is achieved in a finite number of iterations. We next considerthe case when the sequence m p q , m p q , . . . , m p i q , . . . is infinite. The rest of the proof is divided into three partsfor clarity.Convergence. The outline of the convergence proof is as follows. We first demonstrate that the distancebetween any x P S ě| s | X S $ and m p i q or a p i q decreases with every iteration. Using this, we then show that thesequence m p q , m p q , . . . , m p i q , . . . is bounded, and thus, by the Bolzano-Weierstrass theorem and closednessof S ě| s | and S $ , has a subsequence that converges to some m : P S ě| s | X S $ . Referring to the first result again(that the distance between any x P S ě| s | X S $ and m p i q decreases with every iteration), we finally deduce thatthe original sequence m p q , m p q , . . . , m p i q , . . . converges to the same m : P S ě| s | X S $ as any of its infinitesubsequences. S ě| s | X S $ is nonempty by Proposition E.3 . Let us consider some x P S ě| s | X S $ together with m p i q and a p i q taken from the sequences m p q , m p q , . . . , m p i q , . . . and a p q , a p q , . . . , a p i q , . . . for some i ě . Then, we have } x ´ m p i q } “ } x ´ a p i ` q ` a p i ` q ´ m p i q } “ } x ´ a p i ` q } loooooomoooooon ě ` } a p i ` q ´ m p i q } loooooooomoooooooon ě ` ¨ x m p i q ´ a p i ` q , a p i ` q ´ x y loooooooooooooooomoooooooooooooooon ě . (77)The first two terms in the second line of (77) are nonnegative by the definition of the norm. The nonnegativityof the last term follows from the containing-half-space inequality (65) and a p i ` q “ P S $ r m p i q s . Therefore, (77)implies that } x ´ m p i q } ě } x ´ a p i ` q } ` } a p i ` q ´ m p i q } @ i ě (78) For the foundations of AP algorithms in a more general context of arbitrary closed convex subsets of Hilbert spaces, we refer aninterested reader to the seminal works by Bregman [4], Gurin et al. [5], and Boyle & Dykstra [6]. We remind the reader that a p i q “ P S $ r m p i ´ q s for any i ą in the case of the AP-B algorithm. P ALGORITHMS 23 and } x ´ m p i q } ě } x ´ a p i ` q } @ i ě . (79)Replacing m p i q by a p i ` q and a p i ` q by m p i ` q in (77), and using the same argumentation as above, including m p i ` q “ P S ě| s | r a p i ` q s , we deduce that } x ´ a p i ` q } ě } x ´ m p i ` q } ` } m p i ` q ´ a p i ` q } @ i ě ´ (80)and } x ´ a p i ` q } ě } x ´ m p i ` q } @ i ě ´ . (81)The validity of (80) and (81) for not only i ě but also i “ ´ follows from the particular initial conditionsof the AP-B algorithm. Combining (79) and (81) yields } x ´ m p q } ě } x ´ m p q } ě . . . ě } x ´ m p i q } ě . . . (82)and, equivalently, } x ´ m p q } ě } x ´ m p q } ě . . . ě } x ´ m p i q } ě . . . . (83)(83) states that the sequence } x ´ m p q } , } x ´ m p q } , . . . } x ´ m p i q } , . . . is monotonically decreasing. Thissequence is bounded from below by 0 because of the nonnegativity of the norm, and therefore, it converges toits greatest lower bound L ě by the monotone convergence theorem (see Proposition D.10 ). Thus, for every (cid:15) ą , there exists an N p (cid:15) q such that ď } x ´ m p i q } ´ L ď (cid:15) whenever i ą N p (cid:15) q . It follows then from (79)and (81) that L ď } x ´ a p i ` q } ď } x ´ m p i q } ď L ` (cid:15) , so that ď } x ´ m p i q } ´ } x ´ a p i ` q } ă (cid:15) , andbecause of (78), also ď } a p i ` q ´ m p i q } ă ? (cid:15) whenever i ą N p (cid:15) q , i.e., the sequence } a p q ´ m p q } , } a p q ´ m p q } , . . . , } a p i q ´ m p i ´ q } , . . . converges to 0. If the sequence converges, then any of its infinite subsequences(see Definition D.11 ) converges as well, because the removal of elements from the sequence does not changethe validity of the convergence condition: @ (cid:15) ą D N p (cid:15) q : i ą N p (cid:15) q ùñ } a p k i ` q ´ m p k i q } ă (cid:15), (84)where k i ą k j for i ą j ě and k i ě i for i ě .Next, we have from the triangle inequality (63) that } m p i q ´ m p j q } ď } x ´ m p i q } ` } x ´ m p j q } @ i, j ě . (85)Moreover, according to (82), } x ´ m p i q } ď } x ´ m p q } for i ě . Thus, for all i, j ě , } m p i q ´ m p j q } ď ¨ } x ´ m p q } , (86) i.e., the sequence m p q , m p q , . . . , m p i q , . . . is bounded (see Definition D.7 ). Consequently, according to theBolzano-Weierstrass theorem (see Proposition D.12 ), this sequence has a subsequence m p k q , m p k q , . . . , m p k i q , . . . that converges to some m : P R n : @ (cid:15) ą D N p (cid:15) q : i ą N p (cid:15) q ùñ } m : ´ m p k i q } ă (cid:15). (87)We show now that m : P S ě| s | X S $ . By construction, m p i q P S ě| s | for every i ě . According to (87),there exists an m p i q for any (cid:15) ą such that } m p i q ´ m : } ă (cid:15) . Hence, m : is a contact point of S ě| s | (see Definition D.4 ). Moreover, because the latter set is closed (see Definition D.5 and Proposition E.1 ), m : P S ě| s | .Next, by exploiting the triangle inequality (63), we can write } m : ´ a p k i ` q } ď } a p k i ` q ´ m p k i q } ` } m : ´ m p k i q } @ i ě . (88)Combining (88) with (84) and (87) and introducing N p (cid:15) q “ max r N p (cid:15) { q , N p (cid:15) { qs , we get @ (cid:15) ą D N p (cid:15) q : i ą N p (cid:15) q ùñ } m : ´ a p k i ` q } ă (cid:15), (89)i.e., the subsequence a p k ` q , a p k ` q , . . . , a p k i ` q , . . . converges to m : . The set S $ is closed (see Proposi-tion E.2 ), and a p i q P S $ for every i ě by construction. Therefore, using the same argumentation as for thesubsequence m p k q , m p k q , . . . , m p k i q , . . . , we conclude that m : P S $ .Finally, because m : P S ě| s | X S $ , (82) gives } m : ´ m p q } ě } m : ´ m p q } ě . . . ě } m : ´ m p i q } ě . . . . (90)In the light of (90), the statement of (87) generalizes to @ (cid:15) ą D N p (cid:15) q : i ą N p (cid:15) q ùñ } m : ´ m p i q } ă (cid:15), (91)where N p (cid:15) q “ k N p (cid:15) q` . Thus, the sequence m p q , m p q , . . . , m p i q , . . . converges to m : P S ě| s | X S $ .Monotonicity. The monotonicity of the convergence of the sequence m p q , m p q , . . . , m p i q , . . . to m : is declaredby (90).Rate. The key point in establishing the geometric convergence of the AP-B algorithm is the fact that theintersection of S $ and the interior of S ě| s | is nonempty. Using this fact, we first show that the distances between m p i q and S ˝ ě| s | X S $ or a p i ` q and S ˝ ě| s | X S $ can be bounded by, respectively, } m p i q ´ a p i ` q } or } m p i ` q ´ a p i ` q } multiplied by a universal factor that is greater than one and independent of the iteration number i . In the next step,we exploit properties of the metric projection to obtain two additional inequalities, which, in combination with thefirst result, allow deriving a decreasing geometric sequence that bounds } m p q ´ m : } , } m p q ´ m : } , . . . , } m p i q ´ m : } , . . . from above. P ALGORITHMS 25 To start, let us consider some x P S ˝ ě| s | X S $ . Such an element exists according to Proposition E.3 . Also, thereexists an (cid:15) ą such that y P S ě| s | if } x ´ y } ă (cid:15) because x belongs to the interior of S ě| s | (see Definition D.2 ).If so, then it is also possible to choose a positive β ă (cid:15) such that y P S ě| s | if } x ´ y } ď β . We now introduce z p i q “ α i α i ` β ¨ x ` βα i ` β ¨ a p i q , i ě , (92)where α i “ } a p i q ´ P S ě| s | r a p i q s} “ } a p i q ´ m p i q } . Note that α i {p α i ` β q P p , q , and β {p α i ` β q “ ´ α i {p α i ` β q .Moreover, x P S $ , and a p i q P S $ by construction. Therefore, by the Definition D.1 of a convex set and the factthat S $ is convex (see Proposition E.2 ), we have z p i q P S $ . On the other hand, (92) can be rewritten as z p i q “ α i α i ` β ¨ ´ x ` βα i ¨ p a p i q ´ m p i q q ¯loooooooooooooomoooooooooooooon y ` βα i ` β ¨ m p i q . (93)In the above expression, } x ´ y } “ β . Therefore, y P S ě| s | by the definition of β . Moreover, m p i q P S ě| s | by construction, which implies z p i q P S ě| s | , because S ě| s | is convex (see Proposition E.1 ). Hence, altogether, weconclude that z p i q P S ě| s | X S $ for i ě .Based on the above consideration, we show now that } m p i q ´ P S ě| s | X S $ r m p i q s} ď } m p i q ´ a p i ` q } ¨ p ` } x ´ m p q } { β q (94)and } a p i ` q ´ P S ě| s | X S $ r a p i ` q s} ď } m p i ` q ´ a p i ` q } ¨ p ` } x ´ m p q } { β q (95)for i ě . To demonstrate (94), note that } m p i q ´ P S ě| s | X S $ r m p i q s} ď } m p i q ´ z p i q } ď } m p i q ´ a p i ` q } ` } a p i ` q ´ z p i ` q } “ } m p i q ´ a p i ` q } ` α i ` α i ` ` β ¨ } x ´ a p i ` q } ď } m p i q ´ a p i ` q } ` α i ` β ¨ } x ´ m p q } “ } m p i q ´ a p i ` q } ` } m p i ` q ´ a p i ` q } β ¨ } x ´ m p q } ď } m p i q ´ a p i ` q } ` } m p i q ´ a p i ` q } β ¨ } x ´ m p q } . (96)In (96), we used z P S ě| s | X S $ and the Definition D.13 of the projection operator (the first line), the triangleinequality (63) (the second line), (92) (the third line), combined inequalities (79) and (81) (the fourth line), and the Definition D.13 of the projection operator again (the last line). Similarly to (96), we can write } a p i ` q ´ P S ě| s | X S $ r a p i ` q s} ď } a p i ` q ´ z p i ` q } “ α i ` α i ` ` β ¨ } x ´ a p i ` q } ď α i ` β ¨ } x ´ m p q } ď α i ` ` α i ` β ¨ } x ´ m p q } “ } m p i ` q ´ a p i ` q } ` } m p i ` q ´ a p i ` q } β ¨ } x ´ m p q } , (97)which proves (95).Next, we derive two additional inequalities. In particular, we have } m p i q ´ P S ě| s | X S $ r m p i q s} ´ } a p i ` q ´ P S ě| s | X S $ r a p i ` q s} ě } m p i q ´ P S ě| s | X S $ r m p i q s} ´ } a p i ` q ´ P S ě| s | X S $ r m p i q s} “ } m p i q ´ P S ě| s | X S $ r m p i q s} ´ }p a p i ` q ´ m p i q q ` p m p i q ´ P S ě| s | X S $ r m p i q sq} “ ´} a p i ` q ´ m p i q } ` ¨ x a p i ` q ´ m p i q , P S ě| s | X S $ r m p i q s ´ m p i q y“ ´} a p i ` q ´ m p i q } ` ¨ x a p i ` q ´ m p i q , p P S ě| s | X S $ r m p i q s ´ a p i ` q q ` p a p i ` q ´ m p i q qy“ } a p i ` q ´ m p i q } ` ¨ x m p i q ´ a p i ` q , a p i ` q ´ P S ě| s | X S $ r m p i q sy looooooooooooooooooooooooomooooooooooooooooooooooooon ě ě } a p i ` q ´ m p i q } (98)for i ě . Thus, } m p i q ´ P S ě| s | X S $ r m p i q s} ´ } a p i ` q ´ P S ě| s | X S $ r a p i ` q s} ě } m p i q ´ a p i ` q } , i ě . (99)In (98), we used the Definition D.13 of the projection operator (the second line) and the containing-half-spaceinequality (65) along with m p i q “ P S ě| s | r a p i q s (the sixth line). Replacing a p i ` q by m p i ` q and m p i q by a p i ` q and repeating the same steps as in (98), we obtain } a p i ` q ´ P S ě| s | X S $ r a p i ` q s} ´ } m p i ` q ´ P S ě| s | X S $ r m p i ` q s} ě } m p i ` q ´ a p i ` q } , i ě . (100)Combining (94) with (99) and (95) with (100), we obtain, respectively, ˆ ´ p ` } x ´ m p q } { β q ˙ ¨ } m p i q ´ P S ě| s | X S $ r m p i q s} ě } a p i ` q ´ P S ě| s | X S $ r a p i ` q s} (101)and ˆ ´ p ` } x ´ m p q } { β q ˙ ¨ } a p i ` q ´ P S ě| s | X S $ r a p i ` q s} ě } m p i ` q ´ P S ě| s | X S $ r m p i ` q s} , (102) P ALGORITHMS 27 which then lead to ˆ ´ p ` } x ´ m p q } { β q ˙loooooooooooooooooomoooooooooooooooooon r ă ¨} m p i q ´ P S ě| s | X S $ r m p i q s} ě } m p i ` q ´ P S ě| s | X S $ r m p i ` q s} (103)for i ě . Starting with i “ and applying (103) iteratively, we get r i ¨ } m p q ´ P S ě| s | X S $ r m p q s} ě } m p i q ´ P S ě| s | X S $ r m p i q s} , i ě . (104)According to the triangle inequality (63), } m p i q ´ m : } ď } m p i q ´ P S ě| s | X S $ r m p i q s} ` } m : ´ P S ě| s | X S $ r m p i q s} (105)and [see (64)] } m p j q ´ m : } ě ˇˇ } m p j q ´ P S ě| s | X S $ r m p i q s} ´ } m : ´ P S ě| s | X S $ r m p i q s} ˇˇ . (106)(106) and (91) together imply that a sequence } m p q ´ P S ě| s | X S $ r m p i q s} , } m p q ´ P S ě| s | X S $ r m p i q s} , . . . , } m p j q ´ P S ě| s | X S $ r m p i q s} , . . . (107)converges to } m : ´ P S ě| s | X S $ r m p i q s} for every i ě . On the other hand, this sequence is monotonicallydecreasing [see (82)], and thus, by the monotone convergence theorem ( Proposition D.10 ), it converges to itsgreatest lower bound, i.e., } m p j q ´ P S ě| s | X S $ r m p i q s} ě } m : ´ P S ě| s | X S $ r m p i q s} for all i, j ě . Consequently,(105) reduces to } m p i q ´ m : } ď ¨ } m p i q ´ P S ě| s | X S $ r m p i q s} . (108)Applying (108) to (104), we finally obtain γ ¨ r i ě } m p i q ´ m : } , i ě , (109)where γ “ ¨ } m p q ´ P S ě| s | X S $ r m p q s} ą . (cid:4) Remark. 1) The convergence proof of the iterative scheme AP-B relies entirely on the convexity and closednessof the constraint sets. Therefore, this algorithm extends to more general sets than S ě| s | and S $ . 2) The geometricnature of the convergence requires additionally that the interior of at least one of the constraint sets is nonemptyand shares elements with the other set. Proposition III.2. Consider m P M ω and c P C d with | c j | “ ř n { νk “ p ˜ c ν ¨ k ¨ e ı πν p k ´ qp j ´ q{ n q , where ˜ c ν ¨ k P C and n { ν P N . If $ ě ω and ν ě $ ` ω ´ , then a sequence m p q , m p q , . . . , m p i q , . . . formed by the AP-Balgorithm for (cid:15) tol “ and N iter Ñ `8 converges to m .Proof. The proof relies on two auxiliary results that apply to the m and c specified in the proposition: ‚ For every q P R , z “ q ` p| c | ´ q q ˝ θ p| c | ´ q q ùñ z j “ n { ν ÿ k “ ` ˜ z ν ¨ k ¨ e ı πν p k ´ qp j ´ q{ n ˘ , j P I n . (110) ‚ For z defined by (110), P S $ r m ˝ z s “ x z y ¨ m , (111)where, x| z |y “ n ř ni “ | z i | .To show (110), consider a g P R n whose elements form a periodic sequence with the fundamental frequency ν P N such that n { ν P N . Like any element of R n , g can be expressed through its DFT: g j “ ? n n ÿ k “ ` ˜ g k ¨ e ı π p k ´ qp j ´ q{ n ˘ , j P I n . (112)On the other hand, the periodicity of g , g , . . . , g n implies that, for every k P I n , ˜ g k “ p Fg q k “ ? n n ÿ j “ ` g j ¨ e ´ ı π p k ´ qp j ´ q{ n ˘ “ ? n n { ν ÿ j “ g j ν ÿ l “ e ´ ı π pp l ´ q¨p n { ν q`p j ´ qqp k ´ q{ n “ ? n n { ν ÿ j “ ˆ g j ¨ e ´ ı π p k ´ qp j ´ q{ n ν ÿ l “ ` e ´ ı π p k ´ q{ ν ˘ l ´ ˙ “ ˆ ´ e ´ ı πk ´ e ´ ı πk { ν ˙loooooooomoooooooon “ , if p k { ν qR N ¨ ? n n { ν ÿ j “ ` g j ¨ e ´ ı π p k ´ qp j ´ q{ n ˘ . (113)Combining (112) and (113) gives g j “ ? n n { ν ÿ k “ ` ˜ g ν ¨ k ¨ e ı πν p k ´ qp j ´ q{ n ˘ , j P I n . (114)Now, note that | c | defined in the proposition is also periodic with the fundamental frequency ν such that n { ν P N .If so, then the same holds for z defined by (110) because adding a constant or rectifying a function does notchange its periodicity properties. Combining this result with (114) validates the claim of (110). P ALGORITHMS 29 To show (111), consider W $ F p m ˝ z q . For every r P I n , we have p W $ F p m ˝ z qq r “ p W $ q rr ? n n ÿ j “ ` m j ¨ z j ¨ e ´ ı π p r ´ qp j ´ q{ n ˘ “ p W $ q rr ? n n ÿ j “ ˆ m j ¨ „ n ÿ l “ ˆ z l ¨ n n ÿ k “ e ı π p k ´ qp j ´ l q{ n loooooooooooomoooooooooooon δ j,l ˙ ¨ e ´ ı π p r ´ qp j ´ q{ n ˙ “ p W $ q rr ? n n ÿ k “ ˆ ? n n ÿ l “ ` z l ¨ e ´ ı π p k ´ qp l ´ q{ n ˘ ¨ ? n n ÿ j “ ` m j ¨ e ´ ı π p j ´ qp r ´ k q{ n ˘˙ “ p W $ q rr ? n n ÿ k “ ` ˜ z k ¨ ˜ m r ´ k ˘ “ ˜ z ¨ p W $ q rr ? n ¨ ˜ m r ` p W $ q rr ? n n ´ ν ` ÿ k “ ν ` ` ˜ z k ¨ ˜ m r ´ k looooooooooooooomooooooooooooooon “ ðù ω ď $ ď ν ´ ω ` ˘ “ x z y ¨ ˜ m r . (115)When writing the second equality above, we used the orthonormality of F ´ . Combining (115) with the definitionof P S $ r . . . s (see (10) in the main text), we obtain (111).After establishing (110) and (111), consider the sequences m p q , m p q , . . . and a p q , a p q , . . . . Note that m p q “ m ˝ | c |“ m ˝ ` q p q ` p| c | ´ q p q q ˝ θ p| c | ´ q p q q ˘ , (116)where q p q “ . Hence, m p q can be expressed as an elementwise product of the true modulator m and a vectorthat satisfies (110). Let us now assume that, for some i ě , m p i q can be expressed as m p i q “ m ˝ ` q p i q ` p| c | ´ q p i q q ˝ θ p| c | ´ q p i q q ˘ . (117)Then, by (111), a p i ` q “ P S $ r m p i q s “ q p i ` q ¨ m , (118)where, q p i ` q “ x q p i q ` p| c | ´ q p i q q ˝ θ p| c | ´ q p i q qy . (119)Further, by the definition of P S ě| s | r . . . s (see (10) in the main text), m p i ` q “ P S ě| s | r a p i ` q s“ m ˝ ` q p i ` q ` p| c | ´ q p i ` q q ˝ θ p| c | ´ q p i ` q q ˘ , (120)i.e., m p i ` q can also be expressed as the product of m and a vector that satisfies (110). Hence, we conclude byusing mathematical induction that, for every i ě , a p i q “ q p i q m , (121) and q p i q “ x q p i ´ q ` p| c | ´ q p i ´ q q ˝ θ p| c | ´ q p i ´ q qy , (122)with q “ .Next, observe that, by (122), q p i q ´ q p i ´ q “ xp| c | ´ q p i ´ q q ˝ θ p| c | ´ q p i ´ q qyě , (123)i.e., the sequence q p q , q p q , . . . is monotonically increasing. Moreover, it follows from (122) that, for every i ě , q p i ´ q ď ùñ q p i q “ q p i ´ q ` xp| c | ´ q p i ´ q q ˝ θ p| c | ´ q p i ´ q qyď q p i ´ q ` p ´ q p i ´ q q ¨ θ p ´ q p i ´ q q loooooomoooooon “ “ . (124)Taken together with q p q “ , (124) implies that the sequence q p q , q p q , . . . is bounded from above by 1. Hence,by the monotone convergence theorem (see Proposition D.10 ), q p q , q p q , . . . converges to its least upper bound ¯ q ď . If we assume that ¯ q ă , then the convergence of q p q , q p q , . . . to q p q implies that there exists an N such that ¯ q ´ q p N q ď xp| c | ´ ¯ q q ˝ θ p| c | ´ ¯ q qy{ . (125)By (122), q p N ` q ´ q p N q “ xp| c | ´ q p N q q ˝ θ p| c | ´ q p N q qyě xp| c | ´ ¯ q q ˝ θ p| c | ´ ¯ q qyą xp| c | ´ ¯ q q ˝ θ p| c | ´ ¯ q qy{ , (126)which means that q p N ` q ą ¯ q , i.e., the initial assumption that ¯ q ă is incorrect. Therefore, ¯ q “ , and q p q , q p q , . . . converges to 1, i.e., @ (cid:15) ą D N p (cid:15) q : i ą N p (cid:15) q ùñ | ´ q p N q | ă (cid:15). (127)Finally, note that | ´ q p N q | ă (cid:15) ùñ } m } ¨ | ´ q p N q | ă } m } ¨ (cid:15) looomooon (cid:15) ùñ } m ´ m ¨ q p N q } ă (cid:15) ùñ } m ´ a p N q } ă (cid:15) looooooooomooooooooon by (121) , (128)i.e., (127) implies that the sequence a p q , a p q , . . . converges to m . In the light of (81) and (82) in the proof of Proposition III.1 , this result allows concluding that m p q , m p q , . . . also converges to m . (cid:4) P ALGORITHMS 31 AP-A algorithm Proposition III.3. A sequence m p q , m p q , . . . , m p i q , . . . formed by the AP-A algorithm for (cid:15) tol “ and N iter Ñ`8 converges to some m : P S ě| s | X S $ . The convergence is monotonic, i.e., } m p i ` q ´ m : } ď } m p i q ´ m : } for i ě .Proof. Equivalently to the AP-B algorithm, termination of the sequence m p q , m p q , . . . , m p i q , . . . with somefinite i “ N implies that m p N q “ m : P S ě| s | X S $ . Therefore, we next consider the case when the sequence isinfinite. The main idea behind the proof is to show that the inequalities (78) and (80) apply not only to the AP-Bbut also to the AP-A algorithm. When that is established, we can proceed along the path of the convergenceproof of the AP-B scheme.To this end, we first derive some auxiliary (in)equalities. Specifically, it follows from the definition of theoperator P S $ [see (74)] that, for all z , y P R n , x z , P S $ r y sy “ z T F ´ W $ Fy “ z T F ´ W $ W $ Fy “ z T F ´ W $ FF ´ W $ Fy “ p z T F ´ W $ F q ˚ p F ´ W $ Fy q “ p z T FW $ F ´ qp F ´ W $ Fy q“ pp FW $ F ´ q T z q T p F ´ W $ Fy q “ p F ´ W $ Fz q T p F ´ W $ Fy q“ x P S $ r z s , P S $ r y sy . (129)Here, T and * mark, respectively, the transposition and complex conjugation. We used the following propertiesof matrices W $ and F in (129): 1) W $ W $ “ W $ , 2) W ˚ $ “ W T $ “ W $ , 3) F ˚ “ F ´ , 4) F T “ F . Theresult of (129) can be rewritten as x z ´ P S $ r z s , P S $ r y sy “ . (130)(130) is a particular instance of a more general result that the difference between any z P R n and its projectiononto a linear subspace of R n is perpendicular to any element of that subspace. Using (130), we obtain thefollowing: } z } “ x z , z y “ x P S $ r z s ` p z ´ P S $ r z sq , P S $ r z s ` p z ´ P S $ r z sqy“ } P S $ r z s} ` } z ´ P S $ r z s} ě } P S $ r z s} . (131)Applying (131) to the line 6 of the AP-A algorithm ` λ “ } m p i ´ q ´ a p i ´ q } {} b p i q } ˘ , we get λ ě , (132)with the equality holding if and only if p m p i ´ q ´ a p i ´ q q is mapped by P S $ to itself, i.e., m p i ´ q P S $ .However, this would mean that the convergence was reached at iteration i ´ . Finally, we note that line 7 ofthe AP-A algorithm ` a p i q “ a p i ´ q ` λ ¨ b p i q ˘ implies p P S $ r m p i q s ´ a p i ` q q “ p λ ´ q ¨ p a p i q ´ P S $ r m p i q sq (133) and x m p i q ´ a p i ` q , m p i q ´ a p i q y “ x m p i q ´ a p i q , m p i q ´ a p i q y ´ λ ¨ x P S $ r m p i q s ´ a p i ` q , m p i q ´ a p i q qy“ } m p i q ´ a p i q } ´ } m p i q ´ a p i q } } P S $ r m p i q s ´ a p i q } ¨ } P S $ r m p i q s ´ a p i q } “ . (134)In (134) we applied (129) to the term x P S $ r m p i q s ´ a p i ` q , m p i q ´ a p i q qy .We are now ready to prove that (78) and (80) hold for the AP-A algorithm. In particular, we have } x ´ m p i q } “ } x ´ a p i ` q ` a p i ` q ´ m p i q } “ } x ´ a p i ` q } ` } a p i ` q ´ m p i q } ` ¨ x m p i q ´ a p i ` q , a p i ` q ´ x y . (135)The first two terms in the second line of (135) are nonnegative by the definition of the norm. The last term isnonnegative too (note that x P S ě| s | X S $ ): x m p i q ´ a p i ` q , a p i ` q ´ x y “ x m p i q ´ P S $ r m p i q s ` P S $ r m p i q s ´ a p i ` q , a p i ` q ´ x y“ x m p i q ´ P S $ r m p i q s , a p i ` q ´ x y looooooooooooooooooomooooooooooooooooooon “ by (130) `x P S $ r m p i q s ´ a p i ` q , a p i ` q ´ x y“ p λ ´ q ¨ x a p i q ´ P S $ r m p i q s , a p i ` q ´ x y looooooooooooooooooooooooomooooooooooooooooooooooooon by (133) “ p λ ´ q ¨ x a p i q ´ m p i q ` m p i q ´ P S $ r m p i q s , a p i ` q ´ x y“ p λ ´ q ¨ x m p i q ´ P S $ r m p i q s , a p i ` q ´ x y looooooooooooooooooomooooooooooooooooooon “ by (130) `p λ ´ q ¨ x a p i q ´ m p i q , a p i ` q ´ x y“ p λ ´ q ¨ x a p i q ´ m p i q , a p i ` q ´ m p i q ` m p i q ´ x y“ p λ ´ q ¨ x a p i q ´ m p i q , a p i ` q ´ m p i q y loooooooooooooooomoooooooooooooooon “ by (134) ` p λ ´ q loomoon ě by (132) ¨ x a p i q ´ m p i q , m p i q ´ x y looooooooooooomooooooooooooon ě by (65) ě . (136)Therefore, } x ´ m p i q } ě } x ´ a p i ` q } ` } a p i ` q ´ m p i q } . (137)The derivation of the inequality (80) for the AP-A algorithm is equivalent to that for the AP-B: } x ´ a p i ` q } “ } x ´ m p i ` q ` m p i ` q ´ a p i ` q } “ } x ´ m p i ` q } ` } m p i ` q ´ a p i ` q } ` ¨ x a p i ` q ´ m p i ` q , m p i ` q ´ x y loooooooooooooooooomoooooooooooooooooon ě by (65) . (138)Thus, } x ´ a p i ` q } ě } x ´ m p i ` q } ` } m p i ` q ´ a p i ` q } . (139) P ALGORITHMS 33 The rest of the proof follows step by step the proof of Proposition III.1 . Indeed, starting with the inequalities(78) and (80), to which (137) and (139) are equivalent, the proof of Proposition III.1 proceeds based on themand the closedness and convexity of the sets S ě| s | and S $ entirely. The monotonicity of the convergence followsfrom an equivalent of (90), which is derived as a part of the proof of the convergence itself. (cid:4) Remark. 1) Concerning the set S ě| s | , the proof above relies entirely on its closedness and convexity. Thus, theAP-A algorithm is still valid under these more general assumptions. 2) The expressions (129) – (134) applyto projection operators onto any linear space. Hence, the AP-A algorithm still works if S $ is replaced by anarbitrary linear space. Proposition III.4. Consider m P M ω and c P C d with | c j | “ ř n { νk “ p ˜ c ν ¨ k ¨ e ı πν p k ´ qp j ´ q{ n q , where ˜ c ν ¨ k P C and n { ν P N . If $ ě ω and ν ě $ ` ω ´ , then a sequence m p q , m p q , . . . , m p i q , . . . formed by the AP-Aalgorithm for (cid:15) tol “ and N iter Ñ `8 converges to m .Proof. The proof of this proposition goes along the lines of that of Proposition III.2 . First, by using (110)and (111), we derive the result analogous to (121) and (122). Then, we show the monotonic convergence of q p q , q p q , . . . to 1, which assures the convergence of sequences a p q , a p q , . . . and m p q , m p q , . . . to m .To derive the result analogous to (121) and (122), assume that, for some i ě and q p i q P R , a p i q “ q p i q ¨ m , m p i q “ m ˝ ` q p i q ` p| c | ´ q p i q q ˝ θ p| c | ´ q p i q q ˘ . (140)First, note that (140) applies when i “ with q p q “ . Next, by the definition of the AP-A algorithm, we have b p i ` q “ P S $ r m p i q ´ a p i q s“ P S $ r m p i q s ´ a p i q “ P S $ r m ˝ ` q p i q ` p| c | ´ q p i q q ˝ θ p| c | ´ q p i q q ˘ s ´ q p i q ¨ m “ xp| c | ´ q p i q q ˝ θ p| c | ´ q p i q qy ¨ m looooooooooooooooooomooooooooooooooooooon by (110) and (111) , (141) λ “ } m p i q ´ a p i q } } b p i ` q } “ }p| c | ´ q p i q q ˝ θ p| c | ´ q p i q q ˝ m } xp| c | ´ q p i q q ˝ θ p| c | ´ q p i q qy ¨ } m } , (142) a p i ` q “ a p i q ` λ ¨ b p i ` q “ ˆ q p i q ` }p| c | ´ q p i q q ˝ θ p| c | ´ q p i q q ˝ m } } m } ˙looooooooooooooooooooooooooomooooooooooooooooooooooooooon q p i ` q ¨ m , (143) and m p i ` q “ P S ě| s | r a p i ` q s , “ m ˝ ` q p i ` q ` p| c | ´ q p i ` q q ˝ θ p| c | ´ q p i ` q q ˘ . (144)Applying the principle of mathematical induction to the above results, we conclude that, for any i ě , a p i q “ q p i q ¨ m , (145) q p i q “ ˆ q p i ´ q ` }p| c | ´ q p i ´ q q ˝ θ p| c | ´ q p i ´ q q ˝ m } } m } ˙ , (146)with q p q “ .By (146), q p i q ´ q p i ´ q “ }p| c | ´ q p i ´ q q ˝ θ p| c | ´ q p i ´ q q ˝ m } } m } ě , (147)i.e., the sequence q p q , q p q , . . . is monotonically increasing. Moreover, it follows from (146) that, for every i ě , q p i ´ q ď ùñ q p i q “ ˆ q p i ´ q ` }p| c | ´ q p i ´ q q ˝ θ p| c | ´ q p i ´ q q ˝ m } } m } ˙ ď ˆ q p i ´ q ` }p ´ q p i ´ q q ¨ m } } m } ˙ “ , (148)Taken together with q p q “ , (148) implies that the sequence q p q , q p q , . . . is bounded from above by 1. Hence,by the monotone convergence theorem (see Proposition D.10 ), q p q , q p q , . . . converges to its least upper bound ¯ q ď . If we assume that ¯ q ă , then the convergence of q p q , q p q , . . . to q p q implies that there exists an N such that ¯ q ´ q p N q ď ¨ }p| c | ´ ¯ q q ˝ θ p| c | ´ ¯ q q ˝ m } } m } . (149)By (146), q p N ` q ´ q p N q “ }p| c | ´ q p N q q ˝ θ p| c | ´ q p N q q ˝ m } } m } ě }p| c | ´ ¯ q q ˝ θ p| c | ´ ¯ q q ˝ m } } m } ą ¨ }p| c | ´ ¯ q q ˝ θ p| c | ´ ¯ q q ˝ m } } m } , (150)which means that q p N ` q ą ¯ q , i.e., the initial assumption that ¯ q ă is incorrect. Therefore, ¯ q “ , and q p q , q p q , . . . converges to 1. This, as shown in the last paragraph of the proof of Proposition III.2 , implies that a p q , a p q , . . . and m p q , m p q , . . . converge to m . (cid:4) P ALGORITHMS 35 AP-P algorithm Proposition III.5. A sequence m p q , m p q , . . . , m p i q , . . . formed by the AP-P algorithm for (cid:15) tol “ and N iter Ñ`8 converges to a unique m : P S ě| s | X S $ such that } m : } ď } x } for every x P S ě| s | X S $ . The convergenceis monotonic, i.e., } m p i ` q ´ m : } ď } m p i q ´ m : } for i ě .Proof. The proof follows as a corollary of a more general theorem proved for a finite number of closed convexsets in a Hilbert space by Boyle and Dykstra [6, Theorem 2 ]. Specifically, for m p q “ and particular constraintsets S ě| s | and S $ , the sequence m p q , m p q , . . . , m p i q , . . . formed by the algorithm formulated there (the so-calledDykstra’s algorithm) converges to a unique m : P S ě| s | X S $ such that } m : } ď } x } for any x P S ě| s | X S $ .For our purposes, it is thus enough to show that the sequence m p q , m p q , . . . , m p i q , . . . generated by the AP-Palgorithm converges to the same m : .First, we observe that Dykstra’s algorithm formally turns into the AP-P (except the difference in the initialconditions) if we set r “ , denote g i , ” a p i ´ q , g i , ” m p i ´ q , I i , ” c p i ´ q , and assign I i , ” there. Notethat, originally, I i , “ g i , ´ p g i , ´ I i ´ , q and g i , “ P S $ r g i ´ , ´ I i ´ , s . However, due to the linearity of P S $ [see (74)], we have g i , “ P S $ r g i ´ , s ` P S $ r I i ´ , s“ P S $ r g i ´ , s ` P S $ r g i ´ , s looooomooooon “ g i ´ , ´ P S $ r g i ´ , s looooomooooon “ g i ´ , ` P S $ r I i ´ , s“ P S $ r g i ´ , s ` P S $ r I i ´ , s . (151)Applying (151) iteratively, we obtain g i , “ P S $ r g i ´ , s ` P S $ r I , s “ P S $ r g i ´ , s because I , “ byconstruction. Therefore, given the constraint sets S ě| s | and S $ , I i , can indeed be canceled from the Dykstra’salgorithm by setting it to 0 without any consequences to its convergence properties.Further, assuming m p q “ c p q “ , we can verify that Dykstra’s algorithm after the first iteration coincideswith the initialized AP-P algorithm, i.e., the latter is shifted forward by one iteration with respect to the former.Hence, the sequence m p q , m p q , . . . , m p i q , . . . generated by the AP-P converges to the same m : P S ě| s | X S $ asthe analogous sequence formed by the Dykstra’s algorithm initialized with m p q “ .The monotonicity of the convergence of Dykstra’s, and thus the AP-P, algorithms follows from the fact that allterms on the right-hand side of [6, (3.2)] are nonnegative and that each subsequent iteration only adds additionalterms without discarding the old ones. (cid:4) G. L OWER B OUND ON THE C ONVERGENCE E RROR Proof of the statement about lower bound on the convergence error made in Sections III-A and III-B of themain paper are presented here. Proposition G.1. In the case of the AP-B and AP-A algorithms, the convergence error } m p i ´ q ´ m : } {? n isbounded from below by the infeasibility error (cid:15) p i q for any i ě . Proof. According to (79), which applies to both the AP-B and AP-A algorithms, } m p i ´ q ´ m : } {? n ě } a p i q ´ m : } {? n @ i ě . (152)Moreover, by the Definition D.13 of the projection operator and the fact that m : P S ě| s | , we have } a p i q ´ m : } {? n ě } a p i q ´ P S ě| s | r a p i q s} {? n “ } a p i q ´ m p i q } {? n “ (cid:15) p i q @ i ě . (153)Thus, } m p i ´ q ´ m : } {? n ě (cid:15) p i q for all i ě . (cid:4) Remark. Using the Definition D.13 of the projection operator and the fact that m : P S $ , we similarly concludethat } m p i q ´ m : } {? n ě } m p i q ´ P S $ r m p i q s} {? n @ i ě . (154)The latter inequality also applies in the context of the AP-P algorithm. IGNAL CLASSES 37 SIGNAL CLASSES H. T YPES OF W IDEBAND C ARRIERS F OUND IN P RACTICE Wideband carriers encountered in practice fall into three main classes: (quasi-)harmonic , (quasi-)random ,and spike-train . Below, we provide examples of real-world signals featuring these carrier types and applica-tions of their demodulation.The need for demodulating signals formed of harmonic carriers is well recognized in acoustic imaging [7],[8]. There, sinusoidal wavepackets are used as probing signals. However, in many situations, the interactionbetween sound and matter is nonlinear. This makes the returning waves, whose time-dependent amplitude carriesinformation about the imaged object, harmonic. The possibility of efficient and accurate demodulation of signalsof this type would also allow generalizing the probing wave packets themselves from sinusoidal to harmonic.A representative example of quasi-random carriers manifests in surface electromyography. In particular, anelectrical signal detected at the skin surface during the skeletal muscle activity is an amplitude-modulated colorednoise resulting from low-pass filtering of electric pulse trains generated by a large set of conditionally independentmuscle fibers [9, Ch. 5]. The amplitude component of the recorded electrical signal carries information aboutthe force pattern generated by the muscle being studied [10], [11].Probably the most elaborate applications of amplitude demodulation in the context of wideband signals arefound in human speech processing. Speech signals are built of temporarily structured segments of quasi-harmonicand quasi-random carriers [12] that are amplitude-modulated at different timescales [13], [14]. Amplitudedemodulation of broad-band, as well as sub-band speech demodulation, is widely exploited in automatic speechrecognition [15]–[17], hearing restoration [18], [19] tasks, and fundamental studies of the neural mechanisms ofauditory information processing in the brain [20]–[23]. In all these cases, the modulators and carriers convey theinformation about specific aspects of speech, e.g., semantic meaning, associated emotion, or speaker identity,that need to be extracted. Demodulation of signals with mixed quasi-harmonic and quasi-random carrier typesis also known in diagnostic phonocardiography [24], [25]. There, the extracted modulators of heart sounds areused for the detection of events of normal or abnormal functioning of this organ.One of the most popular applications using demodulation of signals built of the spike-sequence type carriersis the Pulse-Code Modulation (PCM) protocol [26]. There, pulses are regularly spaced with specified locationsand have a known constant amplitude and vanishing width. More complex quasi-regular or stochastic pulsesequences of finite width manifest in electric activities of the cardiac muscles and neurons [9, Ch. 6], [27,Ch. 1]. Physiological and diagnostic information contained in these responses is typically associated with themicrostructure and timing of the pulses. However, these recordings often come contaminated by artifacts or otherphysiological signals that slowly modulate the baseline and amplitude of the pulses [28], [9, Ch. 7.1]. Hence, to These properties enable the pulse-coded signals to be demodulated by a simple low-pass filtering procedure [26]. separate the fast cardiac or neural activity (the carrier) from other slowly changing physiological processes orartifacts (the modulator), the raw signal must be demodulated [28]. Finally, as we discuss in Sections IX.A andIX.B of the main text, carriers of the spike-train type also effectively manifest in applications when modulatorrecovery from nonuniformly or sparsely sampled signals of any origin is needed.I. S YNTHETIC T EST S IGNALS The mathematical models and numerical sampling algorithms of synthetic modulators and carriers examinedin this work are described next. Modulators (recovery tests) For signal recovery tests discussed in Section C, modulators were generated by uniformly sampling from M ω .Without loss of generality, we assumed a subset of M ω whose elements’ norms are fixed to ? n . The samplingwas achieved by using a specially adapted version of the Metropolis-Hastings algorithm (see [29, p. 411] for anintroduction to this method). In the concrete, starting with some m p q P M ω , a sequence m p q , m p q , . . . definedby m p i q “ $&% ? n ¨ r p i q {} r p i q } , if r p i q j ě @ j P I n m p i ´ q , otherwise (155)for i ě was generated. Here, r p i q “ m p i ´ q ` F ´ g p i q , (156)and g p i q P C n such that Re ` g p i q ˘ „ N p , σ q , Im ` g p i q ˘ “ , Re ` g p i q j ˘ „ N p , σ q , Im ` g p i q j ˘ „ N p , σ q , ď j ď ω,g p i q j “ ` g p i q n ` ´ j ˘ ˚ , n ` ´ ω ď j ď n,g p i q j “ , ω ` ď j ď n ` ´ ω, (157)and g j KK g k for j ‰ k . We adjusted the parameter σ in (157) to achieve the acceptance rate of the samplingscheme (155) between 0.4 and 0.6. The initial m p q was taken as m p q “ ? n ¨ p g p q ´ min r g p q s ¨ q{} g p q ´ min r g p q s ¨ } . (158)The iterative process generating the sequence m p q , m p q , . . . was terminated upon the first instance of adherenceto the following equilibration criterion of the underlying Markov chain: ω ´ ¨ ω ÿ j “ ˇˇˇˇˇ i ¨ i ÿ l “ p Fm p l q q j ˇˇˇˇˇ ă . ¨ n ω ´ . (159)The corresponding m p i q was then chosen as a sample point from M ω . IGNAL CLASSES 39 Modulators (performance tests) For the performance, convergence, and robustness tests of the demodulation algorithm discussed in Sec-tions IV – VI of the main text, two types of modulators were considered: ‚ Nonstationary Gaussian modulators were produced by transforming a delta-correlated Gaussian processwith a time-dependent low-pass filter. Specifically, m was calculated as m “ w ˝ p m ´ min r m s ¨ q{ max r m ´ min r m s ¨ s , (160)with w being a window “function” (see (169)), and m defined by m i “ n ÿ j “ ´ g p q j ¨ ` F ´ h p i q ˘ i ´ j ¯ , i P I n . (161)In (161), for every i P I n , Re ` h p i q ˘ „ ` P S ω r g p q s ˘ i , Im ` h p i q ˘ “ , Re ` h p i q j ˘ “ ` P S ω r g p j ´ q s ˘ i , Im ` h p i q j ˘ “ ` P S ω r g p j ´ q s ˘ i , ď j ď ω,h p i q j “ ` h p i q n ` ´ j ˘ ˚ , n ` ´ ω ď j ď n,h p i q j “ , ω ` ď j ď n ` ´ ω, (162)with g p j q , j P I ω ´ , being independent samples from the standard n -dimensional Gaussian distribution,i.e., g p j q l „ N p , q for l P I n , and g p j q l KK g p j q k for l ‰ k . ‚ Maximally-uniformly distributed modulators were created by exploiting the NORTA algorithm [30]. Specif-ically, the modulators were calculated as m “ w ˝ p m ´ min r m s ¨ q{ max r m ´ min r m s ¨ s , (163)with w being a window “function” (see (169)), and m given by m “ P S ω r Φ p F ´ r qs , (164)where Φ p . . . q is the cumulative distribution function of the standard Gaussian random variable and r isgiven by r j “ $’’’’’’’&’’’’’’’%a p ‹ ¨ g , j “ b p ‹ j { ¨ p g j ´ ` ı ¨ g j ´ q , ď j ď t p n ` q{ u b p ‹p n ` q{ ¨ g n , j “ p n ` q{ p r n ` ´ j q ˚ , p n ` q{ ă j ď n . (165)In (165), g is the standard n -dimensional Gaussian random vector, and p ‹ “ ? n ¨ p Fc ‹ q ˝ θ p Fc ‹ q , (166) with c ‹ “ , and the remaining components of c ‹ implicitly defined by p F ´ p {? n q j “ ż `8´8 ż `8´8 `` Φ p x q ´ . ˘ ¨ ` Φ p y q ´ . ˘ ¨ φ p x, y | c ‹ j q ˘ dxdy, ď j ď n. (167)In (167), φ p x, y | c ‹ j q is the probability density function of a 2-dimensional Gaussian random variable withzero mean, unit variance, and the correlation between its two elements equal to c ‹ j . Elements of p are givenby p j “ $’’’’&’’’’% n {p ¨ p ω ´ qq , ď j ď ωn {p ¨ p ω ´ qq , n ` ´ ω ď j ď n , otherwise , (168)If all elements of Fc ‹ were nonnegative, the m defined by (164) would have a rectangular power spectrumwith cutoff frequency ω , and m i „ U p , q for every i P I n . However, in reality, such random vector doesnot exist. Therefore, (166) is used to replace Fc ‹ by the closest point in C n that is elementwise nonnegative.On the other hand, this modification expands the power spectrum of the resulting Φ p F ´ r q in (164) beyondthe intended one. The latter discrepancy is canceled by mapping Φ p F ´ r q onto M ω in (164).The window “function” w in (160) and (163) is used to scale the modulator to zero smoothly at the boundarieswith no effect on the remaining points: w i “ $’’’’&’’’’% sin ´ π ¨p i ´ q ¨ n trn ¯ , ď i ď n trn , n trn ă i ď n ´ n trn cos ´ π ¨p i ´ n ` n trn q ¨ n trn ¯ , n ´ n trn ă i ď n @ i P I . (169)In (169), n trn is the length of the transition window. We assumed n trn “ ¨ r f s { ω s in the present work.For simulations discussed in Sections IV – VI of the main text, we used f s “ . ω was set to 15 Hz for nonstationary Gaussian modulators and 20 Hz for maximally-uniformly distributed modulators. The cutofffrequency control parameter $ of the AP algorithms was fixed to 30 Hz . Carriers (recovery tests) For signal recovery tests discussed in Section C, carriers were generated by uniformly sampling from a subsetof C d formed by all spike-train carriers with a fixed number of spikes n s : C n s d “ x P C d : ` ř ni “ I t u p x i q “ n s ˘ ^ ` ř ni “ I t u p x i q “ n ´ n s ˘( (170)The sampling was achieved by exploiting a customized version of the Metropolis-Hastings algorithm. In partic-ular, starting with some q p q P Q such that Q “ x P N n ˚ : ` ř n ˚ i “ x i “ n ´ ˘ ^ ` x i ă d @ i P I n ˚ ˘( , (171) IGNAL CLASSES 41 where n ´ “ n s ¨ d ´ n , and n ˚ “ min t n ´ , n s ´ u , a sequence q p q , q p q , . . . defined by q p i q “ $&% q p i ´ q ` z i ¨ p e p l i q ´ e p k i q q , if ` q p i ´ q ` z i ¨ p e p l i q ´ e p k i q q ˘ P Q q p i ´ q , otherwise (172)for i ě was generated. In (172), e p l q is an element of N n ˚ with its l -th component equal to one and theremaining components equal to zero; l i and k i are independent random numbers from a uniform distribution on I n ˚ ; z i is a random number from a uniform distribution on the set of all integer numbers between ´ σ and ` σ ,where σ is a positive integer chosen to achieve the acceptance rate of the sampling scheme (172) between 0.4and 0.6.The iterative process generating the sequence q p q , q p q , . . . was terminated upon the first instance of adherenceto the following equilibration criterion of the underlying Markov chain: n ˚ ¨ n ˚ ÿ j “ ˇˇˇˇˇ i ¨ i ÿ l “ q p i q j ´ n ´ n ˚ ˇˇˇˇˇ ă . ¨ n ´ n ˚ . (173)The corresponding q p i q was then taken as a sample point from Q . Next, an extended vector ¯q P N n s with ¯ q j “ $&% d ´ q j , if ď j ď n ˚ d, n ˚ ` ď j ď n s (174)was defined. Its components were then randomly permuted to produce another vector ˜q P N n s . The latter wasused to create an r P N n s ` with r j “ η ` j ÿ i “ ˜ q i , j P I n s , (175)where η is a random integer number (the same for all j P I n s ) from a uniform distribution on I n , and thesummation is assumed to be modulo I n . Finally, the carrier c was generated by taking the zero element of R n and setting its components whose indexes are defined by the components of r to one. Carriers (performance tests) For the performance, convergence, and robustness tests of the demodulation algorithm discussed in Sec-tions IV – VI of the main text, four types of carriers were considered: ‚ Nonstationary sinusoid , c “ cos ` π p f t ` ψ q ˘ , (176)with f “ 200 Hz , t “ f ´ s ¨ p , , . . . , n ´ q T , (177)and ψ “ p g ´ min r g sq{ max r g ´ min r g ss , (178) where g “ P S $ r u s u i „ U p , q , i P I n , (179)such that u i KK u j whenever i ‰ j . ‚ Nonstationary harmonic , c “ n f ÿ l “ ´ p { q l ´ ¨ cos ` πlf p t ` ψ l q ` η l ˘¯ , (180)with f “ 85 Hz , n f “ t f s {p ¨ f q u , η l „ U p , q , l P I n f , (181)and ψ l “ η l ` . ¨ p g ´ min r g sq{ max r g ´ min r g ss , (182)where g “ P ω r u s with ω “ and u i „ U p , q , i P I n , (183)such that η i KK η j and u i KK u j whenever i ‰ j . ‚ Nonstationary spikes , c “ ´ θ ` ´ ř n s i “ h ˚ e p r i q ˘ ˝ ` ´ ř n s i “ h ˚ e p r i q ˘ , (184)where e p r i q is the unit vector with all but the r i -th of its components equal to zero, h is defined by h i “ $’’’’&’’’’% e ´p i ´ q { , if ď i ď e ´p n ´ i ` q { , if n ´ ď i ď n , otherwise , (185)and r , r , . . . , r n s is a sequence of different elements of I n generated following r i “ r i ´ ` d i ` η, (186)until r i ´ ď n ´ d i with r “ . In (186), η is a random number from a uniform distribution on I n (thesame for all i ); d i is a random number taken independently from a uniform distribution on I z i at eachiteration, where z i “ r f s {p $ q ¨ p . ` . ¨ sin p π t i qq s , i P I n . (187)Note that, differently from the spike-train carriers generated for the purpose of recovery tests, (184) definessequences of finite-width spikes. ‚ White-noise , c i „ U p´ , q , i P I n , (188)with c i KK c j whenever i ‰ j . ERFORMANCE TESTS 43 PERFORMANCE TESTS J. C ONFIGURATIONS OF D EMODULATION A LGORITHMS FOR P ERFORMANCE A NALYSIS The following configurations of the AP and LDC demodulation algorithms were used for the performanceanalysis in Section IV of the main paper. Window splitting Segment lengths n seg “ t , , . . . , u were examined in the case of the AP demodulation approach.In the case of the LDC demodulation approach, the range of segment lengths was more limited, n seg “t , , . . . , u , to make the simulation times feasible. In order to reduce demodulation errors stemming fromthe window decomposition, signals were split into segments with a particular overlap. On top of that, eachsegment was windowed with the Hann function [31]. For each segment length n seg , different overlap spanswere assumed: n olp “ t , , . . . , n seg { u . The AS-based demodulation was performed only with the originalsignal windows using no decomposition. Control parameters of the AP algorithms Overall, only two parameters are associated with the AP demodulation algorithms: 1) the cutoff frequency ω ,and 2) the number of iterations N iter . In all cases, we fixed ω to 40 Hz . A set of different values of N iter wasconsidered, dependent on the particular algorithm: the range from 1 to 600 for AP-B, 1 to 40 for AP-A, and1 to 6000 for AP-P. It was made sure that the maximum iteration numbers in all these sets allow for reachingdemodulation precision that is high enough not to affect the conclusions of the performance analysis. Control parameters of the quadratic programming solvers for the LDC approach In the case of the LDC approach, all elements of the weighting vector w [see (11)] corresponding to frequenciesbelow the threshold ω were assumed to be zero. The remaining elements of w were set to either of t , u .In the case of the OSQP solver, the following control parameters were tuned: ‚ Linear System Solver: { Suite-Sparse-LDL , MKL-Paradiso } ; ‚ Solution Polishing: t false , true u ; ‚ Warm Starting: t false , true u ; ‚ Absolute Tolerance: t ´ , ´ , . . . , ´ u ; ‚ Relative Tolerance = Absolute Tolerance; ‚ Primal Infeasibility Tolerance = Absolute Tolerance; ‚ Dual Infeasibility Tolerance = Absolute Tolerance; ‚ Maximum Iteration Number: ´ ; ‚ Run Time Limit: 0.The Gurobi solver was tried with these settings: ‚ Method: t , , u ” { primal-simplex , dual-simplex , barrier } ; ‚ Optimality Tolerance: t ´ , ´ , . . . , ´ u ; ‚ Feasibility Tolerance = Optimality Tolerance; ‚ Run Time Limit: 0. K. I MPLEMENTATION ON A C OMPUTER All demodulation algorithms considered in the present work were implemented in C and then interfaced withMATLAB (R2018a) for a large-scale management of different instances and data analysis. The C code wascompiled with GCC (v8.3) using no optimization. Evaluation of the discrete Fourier transform, used in theAS and AP approaches, relied on the FFT implementation of the Intel Math Kernel Library 2019 (update 5).The calculations were performed on a Lenovo “ThinkCentre M910t Tower” desktop computer with an IntelCore i7-7700 processor and Linux Ubuntu 16.04 operating system. In order to minimize the influence of otheroperating system processes on the benchmarking results, one of the four CPU cores was dedicated exclusivelyto the execution of the demodulation program. This was achieved by using the “isolplus” option of the kernelscheduler. All computations were done in single-thread mode. DDITIONAL RESULTS 45 ADDITIONAL RESULTS L. E RRORS OF C ARRIER E STIMATES Here, we discuss findings from the performance and robustness analyses of demodulation algorithms (intro-duced in Sections IV and VI of the main text) in terms of carrier estimates. Analogous to the modulator recoveryerror E m , we use E c “ } c ´ ˆc } {} c } next. Performance tests -1 E c -4 -2 -2 -1 E c -4 -2 -2 -1 E c -4 -2 -1 E c -4 -2 -2 -1 E c -4 -2 T c pu ( s ) AP-PAP-BAP-ALDC-GLDC-O AS-LP AS t (s) -1010 0.025 0.05 t (s) t (s) -101Original AP-B AS-LP0 0.025 0.05 t (s) -101 x sinusoidal A white-noiseharmonic B C D spike-train GFE H Fig. 11. Performance evaluation. A – D : Pareto fronts in the p E c , T cpu q plane for different demodulation algorithms applied to the fourdifferent types of test signals when window splitting is used. The green and brown stars mark the results of, respectively, the AS-LPand AS methods. E – H : Examples of the original carriers considered in the present work (black) and their estimates obtained by usingthe AP-B (red) and AS-LP (green) algorithms. High precision of modulator estimates (see Section IV-C in the main text) and boundedness of ˆ c i between ´ and (see Section VII in the main text) predetermine good quality carrier estimates provided by the APapproach. This view is evidenced by the Pareto fronts in the p E c , T cpu q plane shown in Fig. 11 A–D andcomparisons of exemplary c and ˆc in Fig. 11 E–H. We observed noticeable discrepancies between c and ˆc only locally, around points with modulator levels very close to zero (see the signal segment at t “ . inFig. 11 E). That finding is explained by the fact that, in our case, ˆ c i ´ c i “ s i ¨ p ˆ m ´ i ´ m ´ i q . Hence, evensmall differences between m i and ˆ m i on the absolute scale can give notable differences between c i and ˆ c i if m i { ˆ m i " . However, the condition ˆ m i ě s i ( ùñ | c i | ď ) inherent to all ˆc obtained by the AP algorithmsassures tolerable distortions of the carrier estimates even around the points with vanishingly small m i . Remember that c i { ˆ c i “ m i { ˆ m i in our case. The situation is different in the case of the AS-LP approach. Then, | s i { ˆ m i | " , i.e., | ˆ c i | " , are possible(see the signal segment at t “ . in Fig. 11 E). These divergences noticeably increase the recovery errors E c even for sinusoidal signals that AS-LP recovers sufficiently well outside the segments of vanishing m i (seeFig. 11 A, E). Globally-inaccurate recovery of other carrier types demonstrated by the AS-LP (see Fig. 11 B– D,F– H) is predetermined by inaccurate estimates of the modulators (see Section IV-C in the main text). We notethat, for sinusoidal carriers, appropriate estimates ˆc can be obtained by using the original AS method insteadof the AS-LP (Fig. 11 A). However, the AS gives very erroneous modulator estimates ˆm , and hence, does notimprove carrier predictions ˆc considerably, for other types of signals (data not shown). Robustness tests P(0) -2 -1 P(0) -2 -1 P(0) -2 -1 P(0) -2 -1 E c AP-B AP-PAP-AAS-LP sinusoidal A white-noiseharmonic B C D spike-train Fig. 12. Robustness evaluation. A – D : Dependence of the carrier recovery error E c on P p q (the probability of missing points) for thefour types of test signals and different AP algorithms at (cid:15) tol “ ´ (color coding). Fig. 12 shows demodulation results of the AS-LP, AP-B, AP-A, and AP-P algorithms in terms of carrierrecovery error for test signals corrupted by the multiplicative Bernoulli- t , u noise (see Section VI). Theobtained E c vs. P p q relations for the AP algorithms are analogous to their counterparts E m vs. P p q shownin Fig. 5. In contrast, the AS-LP is even more inferior to the AP approach in terms of carrier reconstruction(Fig. 12) than it is in terms of modulator recovery (Fig. 5 A–D). This fact is explained by the presence of the | ˆ c i | " divergences discussed above and illustrated by Fig. 11 E.M. C ONVERGENCE R ATES AT D IFFERENT S IGNAL L ENGTHS The convergence results shown in Fig. 4 of the main paper represent only signals with the length n fixed to sample points. Therefore, we performed additional simulations with different n values to test the impactof this parameter on the rate of the iterative process. We found no clear dependence of N iter required toachieve a particular infeasibility error value (cid:15) on n , except transitional changes due to diminishing contributionsof the boundary effects in some cases (see Fig. 13 C – C ). These findings suggest that the increased T cpu ofdemodulation for longer n is primarily determined by the increased computational demands of single projections(see Section III-D in the main paper). DDITIONAL RESULTS 47 n n n n n n N i t e r AP-B tol =10 -3tol =10 -6tol =10 -9 n n n N i t e r AP-A tol =10 -3tol =10 -6tol =10 -9 n n n N i t e r AP-P tol =10 -2tol =10 -3tol =10 -4 sinusoidal A white-noiseharmonic B C D spike-train C B A D C B A D Fig. 13. Convergence analysis of the AP algorithms at different n and fixed f s “ . A – D : Dependence of the iteration number N iter necessary to reach a specific infeasibility error (cid:15) on the length of the signal for the AP-B algorithm applied to four differentclasses of test signals. Filled circles, rectangles, and triangles label curves for (cid:15) levels of, respectively, ´ , ´ , and ´ . A – D :The same as A – D but shown for the AP-A algorithm. A – D : The same as A – D but shown for the AP-P algorithm and differentlevels of the infeasibility error (cid:15) . t (s) t (s) x t (s) 01 0 0.1 0.2 0.3 t (s) sinusoidal A white-noiseharmonic B C D spike-train Fig. 14. Repeated demodulation of inferred carriers. Gray – sinusoidal ( A ), harmonic ( B ), spike-train ( C ), and white-noise ( D ) carriersinferred by demodulating test signals shown in Fig. 2 A–D with the AP-A algorithm. Black – target repeated modulators following fromthe assumption that the inferred carriers are fully demodulated. Orange – the real repeated modulators obtained with AP-A. N. R EPEATED D EMODULATION Fig. 14 shows typical results of redemodulation of carriers inferred from the four types of synthetic test signalsconsidered in the present work by using the AP-A algorithm. Redemodulation of the carriers returns the identity modulator to an excellent approximation, implying nearly perfect demodulation achieved by the AP approach,as discussed in Section VII of the main paper.O. D EMODULATION OF S PEECH S IGNALS In this section, we present results of additional simulations used to support the statements about the suitabilityof the AP approach to demodulate wideband speech signals in Section VIII of the main text. Synthetic modulators The first question that we considered was up to which values of the cutoff frequency ω the natural speechcarriers meet the recovery condition r n { d s ě ω ´ . As mentioned in the main text, these carriers are ofquasi-random and quasi-harmonic origins, possibly featuring frequency glides. It is well known that typicalfundamental frequencies ( f ) of male and female speaker voice are, respectively, 120 and 210 Hz (see [32] andTable 1 therein). Hence, at least for the harmonic components, the condition r n { d s ě ω ´ is expected to besatisfied with ω ď 60 Hz . That is more than sufficient for appropriate demodulation, assuming that strictly allspectral energy of the speech amplitude modulator is located below 20 Hz (see Section VIII-A). The situationwith the quasi-random and mixed components of speech signals is less certain and must be tested numerically.To this end, we considered four speech-carriers generated by a male speaker uttering prolonged ( „ s)[ α :], [ ş ], and [z], as well as vocable [wa f ] (see Fig. 15 A – A ) at f “ 120 Hz without noticeable amplitudemodulation. The [ α :] is predominantly harmonic, [ ş ] is quasi-random, [z] has a mixed wave-shape, and [wa f ]is quasi-harmonic but with upward and downward frequency glides. These characteristics are further revealedby the power-spectral-density (PSD) plots for the [ α :], [ ş ], and [z] (see Fig. 15 B – B ), and a spectrogram forthe [wa f ] (Fig. 15 B ). We then created test signals as products of mentioned carriers and maximally-uniformlydistributed synthetic modulators (see Section I) with ω “ 25 Hz and f s “ . .Demodulation of the considered test signals with the AP-A algorithm allowed us to achieve high-accuracymodulator recovery, as shown in Fig. 15 C – C (note the insets with E m and E c values there). In more detail,we found that, for the [ α :], [ ş ], [z], and [wa f ], the distances between carrier points with absolute values of,respectively, ě . , ě . , ě . , and ě . were below that needed to satisfy r n { d s ě ω ´ at ω “ 20 Hz . For the [ ş ], [z], and [wa f ], 80 % of the required points were ě . . The demodulation qualityremained reasonably good when ω was increased to 50 Hz , resulting in E m values of ¨ ´ , ¨ ´ , ¨ ´ ,and ¨ ´ for, respectively, [ α :], [ ş ], [z], and [wa f ] carriers. Natural modulators Next, we aimed to clarify whether the AP approach can properly separate m and c of speech signals usingthe dynamic range compression discussed in Section VIII-B of the main text. For this purpose, we considered DDITIONAL RESULTS 49 t (ms) -101 c t (ms) -101 c t (ms) -101 c t (ms) -101 c f (kHz) f (kHz) f (kHz) t (s) f ( k H z ) -80-60-40-20 PSD (dB / Hz) B B B B prolonged [ ʃ ] prolonged [ z ] vocable [ wa ʊ ] prolonged [ ɑ : ] A A A A Original AP-A E m =1 10 -3 E c =1 10 -2 t (s) x E m =3 10 -2 E c =7 10 -2 t (s) x E m =1 10 -2 E c =5 10 -2 t (s) x E m =2 10 -2 E c =7 10 -2 t (s) x C C C C Fig. 15. Demodulation of signals built of synthetic modulators and natural speech carriers. A – A : fragments of carrier-signals ofprolonged steady [ α :] ( A ), [ ş ] ( A ), and [z] ( A ), as well as an extended vocable [wa f ] ( A ). B – B : periodogram estimators of thepower spectral densities of the carriers illustrated in A – A ; B : spectrogram of the vocable [wa f ]. C – C : exemplary segments ofamplitude-modulated carriers from panels A – A , their modulators (black), and modulator estimates obtained with the AP-A algorithm.Insets of panels C – C display values of the modulator and carrier recovery errors. synthetic time series built of the four natural carriers [ α :], [ ş ], [z], and [wa f ] introduced above and modulatorestimate ˆm ˚ of an utterance “. . . protein which forms p . . . ” from Section VIII-B of the main text. We thendemodulated the resulting signals by using the AP-A algorithm combined with the dynamic range compression.The obtained estimates ˆˆm ˚ were in good agreement with the original, as shown in Fig. 16 (note the insets with E m and E c values there). The AP-B, AP-P, and LDC algorithms returned very similar modulator estimates butneeded much longer computing times that mirror the performance analysis presented in Section IV-C of themain paper. Original AP-A E m =9 10 -3 E c =1 10 -2 t (s) m E m =3 10 -2 E c =5 10 -2 t (s) m E m =5 10 -2 E c =6 10 -2 t (s) m E m =5 10 -2 E c =1 10 -1 t (s) m prolonged [ ʃ ] prolonged [ ɑ : ] A B prolonged [ z ] vocable [ wa ʊ ] C D Fig. 16. Demodulation of signals built of natural speech modulators and carriers. A – D compares original modulators ˆm ˚ (black) andtheir estimates ˆˆm ˚ obtained by using the AP-A algorithm combined with the dynamic range compression (orange) for, respectively, [ α :],[ ş ], [z], and [wa f ] carriers. Insets display values of the modulator and carrier recovery errors. R EFERENCES [1] F. Zhang, Matrix Theory . Springer, 2011.[2] A. N. Kolmogorov and S. V. Fomin, Introductory real analysis . Courier Corporation, 1975.[3] E. øCinlar and R. J. Vanderbei, Real and Convex Analysis . Springer, 2013.[4] L. M. Bregman, “The relaxation method of finding the common point of convex sets and its application to the solution of problemsin convex programming,” USSR Comput. Math. & Math. Phys. , vol. 7, no. 3, pp. 200–217, 1967.[5] L. G. Gubin, B. T. Polyak, and E. V. Raik, “The method of projections for finding the common point of convex sets,” USSRComput. Math. & Math. Phys. , vol. 7, no. 6, pp. 1–24, 1967.[6] J. P. Boyle and R. L. Dykstra, “A method for finding projections onto the intersection of convex sets in Hilbert spaces,” in Advancesin Order Restricted Statistical Inference , ser. Lecture Notes in Statistics. Springer, 1986, pp. 28–47.[7] V. F. Humphrey, “Nonlinear propagation in ultrasonic fields: measurements, modelling and harmonic imaging,” Ultrasonics , vol. 38,no. 1, pp. 267–272, 2000.[8] F. A. Duck, “Nonlinear acoustics in diagnostic ultrasound,” Ultrasound Med. Biol. , vol. 28, no. 1, pp. 1–18, 2002.[9] L. Sornmo and P. Laguna, Bioelectrical Signal Processing in Cardiac and Neurological Applications . Academic Press, 2005.[10] G. L. Gottlieb and G. C. Agarwal, “Filtering of electromyographic signals,” Am. J. Phys. Med. Rehab. , vol. 49, no. 2, p. 142, 1970.[11] T. J. Roberts and A. M. Gabald´on, “Interpreting muscle function from EMG: lessons learned from direct measurements of muscleforce,” Integr. Comp. Biol. , vol. 48, no. 2, pp. 312–320, 2008.[12] J. E. Shoup and L. L. Pfeifer, “Acoustic characteristics of speech sounds,” in Contemporary Issues in Experimental Phonetics .Academic Press, 1976, pp. 171–224.[13] R. E. Turner, “Statistical models for natural sounds,” Ph.D. dissertation, University College London, 2010. [Online]. Available:http://discovery.ucl.ac.uk/19231/[14] A. Keitel, J. Gross, and C. Kayser, “Perceptually relevant speech tracking in auditory and motor cortex reflects distinct linguisticfeatures,” PLOS Biol. , vol. 16, no. 3, p. e2004473, 2018. EFERENCES 51 [15] B. E. D. Kingsbury, N. Morgan, and S. Greenberg, “Robust speech recognition using the modulation spectrogram,” Speech Commun. ,vol. 25, no. 1, pp. 117–132, 1998.[16] S. Wu, T. H. Falk, and W.-Y. Chan, “Automatic speech emotion recognition using modulation spectral features,” Speech Commun. ,vol. 53, no. 5, pp. 768–785, 2011.[17] B. Lee and K.-H. Cho, “Brain-inspired speech segmentation for automatic speech recognition using the speech envelope as atemporal reference,” Sci. Rep. , vol. 6, p. 37647, 2016.[18] B. S. Wilson, C. C. Finley, D. T. Lawson, R. D. Wolford, D. K. Eddington, and W. M. Rabinowitz, “Better speech recognitionwith cochlear implants,” Nature , vol. 352, no. 6332, pp. 236–238, 1991.[19] F.-G. Zeng, K. Nie, G. S. Stickney, Y.-Y. Kong, M. Vongphoe, A. Bhargave, C. Wei, and K. Cao, “Speech recognition with amplitudeand frequency modulations,” PNAS , vol. 102, no. 7, pp. 2293–2298, 2005.[20] R. V. Shannon, F.-G. Zeng, V. Kamath, J. Wygonski, and M. Ekelid, “Speech recognition with primarily temporal cues,” Science ,vol. 270, no. 5234, pp. 303–304, 1995.[21] Z. M. Smith, B. Delgutte, and A. J. Oxenham, “Chimaeric sounds reveal dichotomies in auditory perception,” Nature , vol. 416, no.6876, pp. 87–90, 2002.[22] P. X. Joris, C. E. Schreiner, and A. Rees, “Neural processing of amplitude-modulated sounds,” Physiol. Rev. , vol. 84, no. 2, pp.541–577, 2004.[23] U. Goswami, “Speech rhythm and language acquisition: An amplitude modulation phase hierarchy perspective,” Ann. N. Y. Acad.Sci. , vol. 1453, no. 1, pp. 67–78, 2019.[24] A. A. Sarkady, R. R. Clark, and R. Williams, “Computer analysis techniques for phonocardiogram diagnosis,” Comput. Biomed.Res. , vol. 9, no. 4, pp. 349–363, 1976.[25] D. Gill, N. Gavrieli, and N. Intrator, “Detection and identification of heart sounds using homomorphic envelogram and self-organizingprobabilistic model,” in Computers in Cardiology , 2005, pp. 957–960.[26] B. Oliver, J. Pierce, and C. Shannon, “The philosophy of PCM,” Proc. IRE , vol. 36, no. 11, pp. 1324–1331, 1948.[27] F. Rieke, D. Warland, R. de Ruyter van Steveninck, and W. Bialek, Spikes: exploring the neural code . A Bradford Book, 1999.[28] J. Felblinger and C. Boesch, “Amplitude demodulation of the electrocardiogram signal (ECG) for respiration monitoring andcompensation during MR examinations,” Magn. Reson. Med. , vol. 38, no. 1, pp. 129–136, 1997.[29] L. Wasserman, All of statistics: A concise course in statistical inference . Springer, 2004.[30] M. C. Cario and B. L. Nelson, “Modeling and generating random vectors with arbitrary marginal distributions and correlationmatrix,” Department of Industrial Engineering and Management Sciences, Northwestern University, Tech. Rep., 1997.[31] G. Sell and M. Slaney, “Solving demodulation as an optimization problem,”