Okay, I have had some time to think about this, and I think I have solved at least one of the problems I laid out in the original post

>Although I have proposed this algorithm, I am yet to implement it. The main challenge is how the optimization taking into account the heuristics should be performed. The simplest idea would be to optimize each row of the matrix individually using any appropriate technique (e.g. gradient descent, Simultaneous Perturbation Stochastic Approximation). Then repeat this many times. However this would presumably be very inefficient. Besides, we should ideally be optimizing globally or at least using regions larger than a single sample. Another issue is that optimization procedure may change the deltas to minimize the heuristics, instead of distributing the residual in the way that minimizes the heuristics' cost.

The solution to the issue of the heuristics making changes in a way different than the residual dictates, or beyond what they dictate, is as follows. The solution is to enforce the requirement V/S >= C. Where V is the sum of a row, and S is the sum of the absolute values of a row. C is a constant. To enforce the requirement I laid out originally (i.e. that the residual energy must be just distributed, so no kinds of 'accounting' is allowed where values of the opposing sign could also result to minimize heuristics), C would be equal to one. On the other end, to allow some amount of 'creativity', C could be set to be some value 0<C<=1.

In an optimization based approach, these hard limits may be a problem. Thus I propose using the following differentiable cost function: c_2*e^(-v*c_1). Where v is a value that becomes negative when the condition is not met, and c_1 and c_2 are weighting constants. For the proposed requirement, v=V/S-C, or v=(V/S-C)*S if we want it go be scaled also by the magnitude of the changes.

For the optimization algorithm, I am still not completely sure, but at least at first, I think I will go with a global gradient descent. Of course, this will mean it would probably be non-real-time. However, for the main application I am working on, a single voice synthesizer, this could be fine if the analysis is done during the database creation, since the gradient descent doesn't have to be run for synthesis.

I have also thought of another approach, actually that doesn't use an optimizer at all. This approach probably wouldn't work as well as the optimizer-based approach, but it may be able to run I real-time. The idea is to center a window function on each of the voice pulse onsets. Then, the RVPM residual is normalized according to the sum of all of those window functions. Finally, the residual is added to each voice pulse, according to its corresponding window.

A parametric window like the Kaiser-Bessel makes sense to find the best parameter, which should presumably be constant since the window size is fixed. Besides, the performance should not be much of a window since we can just construct a large window once and compute all others using cubic spline interpolation of that window. I have already tested this approximation and shown that it yields an error of less than -200dB for a 2049-point Kaiser-Bessel 2.5 window. Additionally, this window can be modified based on the fact that a real voice pulse should be zero before the start of the opening phase. We can determine a safe maximum opening phase duration, and then set the window to ramp to down to zero before that start.