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.

The #VectorAutoRegression stuff I've been doing can be summarized as
$$ y_t = \sum_{i=1}^p A_i y_{t - i} $$
where each $y_t$ is a $D$-vector and each $A_i$ is a $D \times D$ matrix.

Given an input series of $y$ values, the $A_i$ can be calculated by #LeastSquares minimization:

$$
Y = [ y_p ; y_1; \ldots ; y_{T - 1} ] \\
X = [ y_{p-1}, y_{p-2}, \ldots, y_0 ; y_p, y_{p-1}, \ldots, y_1 ; \ldots ; y_{T-2}, y_{T-3}, \ldots, y_{T-1 - p} ] \\
A = \left(X^T X\right)^{-1} X^T Y
$$
where the matrices have these initial dimensions:
Y : (T - p) × D
X : (T - p) × (p × D)
A : (p × D) × D
then reshape $A$ to $p × (D × D)$

In my case all values are complex (and ^T is conjugate transpose) because each $y$ is an FFT of a block of input data - I'm using FFT size $256$ (making $D = 129 = 256/2+1$ unique bins for real input) overlapped 4x.

More research on #VectorAutoRegression got me digging into #SpectralRadius which is the magnitude of the largest-magnitude #Eigenvalue. This is analogous to the pole radius in regular single variable #ZTransform representation: if it's less than 1 all should be fine, bigger than 1 and it becomes unstable.

So I'm now normalizing all the feedback coefficient matrices by the largest spectral radius among them and a bit more, so that the new largest radius is less than 1, and it seems stable even in the presence of morphing.

The attached is heavily dynamics compressed, as it was a bit peaky otherwise.

I think I got interpolation between feedback coefficient matrices working (thanks to SciPy, which has polar decomposition, log and exp for complex matrices) but it blew up almost immediately, so I'm re-adding the normalization code in the hopes of getting some output that's more than a brief click followed by NaNs.

Trying to make morphing drones modulated by a third signal or something, not sure what it'll turn into...

#VectorAutoRegression