Went back to the #HaarWavelets stuff and found a way to control the rhythmicality of the noise by scaling the variances of the seeding decorrelated gaussian random numbers. With sigma = 1 it's rhythmic, smaller is more uniform. Maybe I could control it per octave instead of all octaves together in unison.

Now I'm thinking of combining this with the FFT-based stuff: X would be the 10D (real; usually 11D but I think it's maybe best to exclude DC for this) energy per octave analysis of the same audio that Y is the 128D (complex, excluding DC) FFT analysis of, then I can simulate X with #VectorAutoRegression like the attached (without generating audio) and feed that into the #VARMA to get Y for audio output.

First successful #VARMA test with small P, Q. (Follows about 25 unsuccessful tests.)

Audio attached has 3 parts, the first is the target sound (Y), the second is the input sound (X), the third is the VARMA model applied to X (normalized in volume) which should hopefully sound like Y.

P = 0 (one matrix for input mapping, no real moving average here...)
Q = 1 (one matrix for auto-regressive feedback)
D = 129 (number of dimensions for FFT vectors)
T = 18103 (number of input vectors)

Computed #SpectralRadius of the feedback matrix was 0.9144670749524257, less than 1 so it's stable (good).

The peak level of the output audio was 0.148676, significantly quieter than I expected (most of my previous tests using Python statsmodels VAR had gains around 1.0e+7 or so...).

#VectorAutoRegression is a pure feedback model, similar to poles in the pole-zero representation of z-transform filters. To add zeros to the model is quite simple:

$$ y_t = \sum_{i=1}^p A_i y_{t - i} + \sum_{i=0}^q B_i x_{t-i}$$

and then add the $x$ vectors and $B$ matrices to the least squares stuff. I think the jargon for this is #VARMA, for Vector Auto-Regressive Moving Average.

I worked it through and coded up most of it before realizing a #FatalFlaw : I need the $x$ input data corresponding to the $y$ output data to do the regression to estimate the $A$ and $B$ parameter matrices, while so far I've just been using WAV files as $y$ series, which means I don't have any $x$.