Merging pattern for FFT solutions I have the generic problem of analysing a discretly sampled signal, lets
say at 100Hz, and would like to compute features (ML application) of this signal at different scales:

An instantiation of this problem is computing the FFT on real-signal with
N samples:

Original signal: X[1],...,X[N]

Suppose I have already computed the FFT of the following
subarrays:

X[N-k+1],...,X[N]

X[N-2*k+1],...,X[N-k]

X[N-3*k+1],...,X[N-2*k]

X[N-4*k+1],...,X[N-3*k]

which is each of size k. Then I would like to compute the FFT
of the subarray

X[N-4*k+1],...,X[N] 

(which is of size 4*k), WITHOUT reconstructing the original array X[.] and
running an ordinary FFT on the array of size 4k. Roughly I would like to do O(N) primitive operations in the best case if that is at all possible.

Is there some recursive property I could use?

PS: Note that if the goal of computing the FFT is substituted with e.g. finding the mean or the variance there exist easy recursive properties that allow merging of the results in constant time.
You can compute the DFT of the total signal from the DFTs of its parts but this is not efficient. What you would need to do is interpolate each DFT to the length of the total array (which is equivalent to zero-padding the sub-blocks to length $N$ before computing the individual FFTs), and then combine the interpolated DFTs with appropriate complex factors. This is computationally less efficient than going back to the time domain and computing the FFT of the whole sequence.

There is no recursion involved (unlike your examples of recursive mean and variance computation), simply because the data blocks are not related to each other. If you were to have a sliding window with a fixed length and you wanted to recompute the DFT after shifting the window by one sample, then you would indeed get a recursive formula for the DFT. This is known as the sliding DFT.