Update jitter noise models

We currently use white jitter noise (in angle), for MOSA and SC. The allocation is not white, and have some relaxation at low frequencies, with a noise in

$$A \Big[1 + \Big( \frac{0.8 \text{mHz}}{f} \Big)^4 \Big]$$

This corresponds to feature (7) in the list, with high priority, deadline January 2023. Note that the noise shape will probably need to be updated in the near future, so that we should allow for pole-zero model parametrization.

This is implemented for LISANode (one of the many changes in https://gitlab.in2p3.fr/j2b.bayle/LISANode/-/merge_requests/350, to be split into its own self-contained merge request). We use a Franklin noise generator. We don’t have Franklin generators in LISA Instrument, so we can either implement it (a bit more work), or use ad-hoc models for now.

changed milestone to %v1.1

added feature label

changed the description

(Nevermind this comment, CBE = allocation right now for some reason)

Is allocation the correct choice, or do we use CBE? Obviously allocation is conservative but it's also often unphysical in spectral character. Right now we have the number from the L3 proposal (10^-8 [m/rt(Hz)]) which was really neither of those specific concepts, but probably closer to the CBE.

Silly question: do we have a CBE for this?

Usually (for now) we implement the CBE (or a simplified version of it), and fall back to allocation if there's no CBE.

Yeah, we don't in the perf model itself, and I think upon a few seconds contemplation that's because we have two different SC designs from the Primes, so we actually can't have one CBE. My bad. Anyhow, the flat jitter level is lower by a factor of 2 than the old value, [5E-9.*sqrt(1+(0.8E-3./f).^4)], and the low frequency ramp up happens as TTL gets less relevant, so as above I rescind my comment!

changed milestone to %v1.2

changed the description

@mastaa, since you implemented the Franklin generator in LISANode, how difficult do you think it is to port it to Python?

It's not gonna be easy but fortunately I have a somewhat working implementation in Python! (As I first did a quick and dirty implementation in Python to then ported it to C++)

One difficulty I see is that in the propagation step (going from time step

$i$

to

$i+1$

) we need to form the matrix-vector product of a constant matrix and the last intermediary noise vector. This translates really well in LISANode since we simulate sample-by-sample but for LISA Instrument where we simulate the full time series I don't know how to efficiently port that to Numpy syntax?

Indeed. Is there a way to define the entire series of intermediary noise vectors, and then use convolution?

Hm, not sure. That's basically the two core equations:

$\vec{\xi}(t_k)$

is the intermediary noise vector which needs to be propagated every time step where

$\vec{w}$

is a vector of Gaussian independent noise (different for different

$k$

). The involved matrices are constant.

It seems not then… Each

$\xi(t_k​)$

depends on the previous one. So that looks like a loop…

So we might go for the easy solution: using ad-hoc composed powerlaws!

Yes, I would also vote for that! Another option could be to start with a continuous time pole zero model and use the Bilinear Transform to design an approximate IIR filter for white noise. Have you tried that before? It seems quite easy and possibly even equivalent to what you propose?

Yes, that's pretty much what we do… except that IIR filters are not super well-behaved. Once has to be careful to implement them as order-2 cells, and they are a pain because we cannot (in theory) chunk the simulation in time anymore.

I mean, we can, but each chunk depends now on the previous one, so everything has to be sequential.

Isn't that true already since we do integration?

You mean, for noise generation?

Yes, that's when we have models with poles. Like 1/f noises.

I also mean for the classical example of calculating the timer deviations by integrating the clock noise. So there are already elements in the simulation which prevent chunking and parallel execution, right?

Probably, yes.

Parallel execution is one thing… not being able to chunk is another!

Ah, really? Operating the IIR with order 2 cells doesn't allow chunking anymore?

Does it? You need the result of the entire previous chunk, right?

In my understanding to calculate the next sample for an IIR filter (i.e. y[n]) you need some x[n], x[n-1], ..., x[n-M] and y[n-1], ..., z[n-N] with a finite M and N depending on the order of the filter. So you need those samples from the previous chunk not the entire one. Where is my loop hole?

You're right, theoretically.

I think I'm projecting a Dask-like implementation of this. And in this case, you map computations on overlapping-chunked arrays, such that you have to be able to express your result as y[n:n+N] = f(x[n:n+M]), where M is the chunk length, including margins. With just this framework, you cannot, right?

I know that Dask has some algorithms to implement those kinds of computations (and even partially parallelize them), but I'll need to take a closer look at them. See https://docs.dask.org/en/stable/generated/dask.array.cumsum.html.

I see your point...

assigned to @j2b.bayle

created merge request !136 (merged) to address this issue

mentioned in merge request !136 (merged)

mentioned in commit 1eed3564

closed with merge request !136 (merged)