How do I determine stationarity from a set of 50 complex values collected every 10 minutes?

Question

I am trying to determine stationarity from a somewhat stochastic process. Every 10 minutes, I collect a set of 50 FFTs, i.e., 1 trial over $50$ seconds, so an FFT occurs every time second. I understand the definition of weak-stationarity with the covariance only being dependent upon the lag interval, the mean and variance remains constant. At the end of the day I want to answer the question "when do I stop collecting data?" i.e. my data is "weakly-stationary". Quick preface, I have been working on this problem for weeks now and am at somewhat of a loss.

Further information on my data, and known parameters:

• Each FFT only outputs one complex value which corresponds to the amplitude and phase of a desired frequency signal.
• Sample rate for 1 FFT is known (sample rate of x(t))
• Bandwidth is known for the FFT calculation
• Sampling interval known of $x(t)$
• Spectral lines is known
• Frame size is known
• Frequency resolution known
• Block size known
• Frequency of original signal is known
• Window type is Hann with linear spectral format
• I know what the time signal should look like.
• Assume data comes from a signal that is i.i.d.
• I only have 1 complex value for an FFT since the computation is only interested in outputing that value from that desired frequency.

What I am unsure about my data:

• I’m 95% confident that my real and I’m amplitudes for the complex values went through an RMS calculation. Maybe ignore this point, just adding info for more context.

What is unknown:

• The data of the original signal is not known.
• If my complex values are for a one or two spectral lines. I noticed that my data are split into two regimes, separated by some small distance/phase. I am unsure how to account for this.

My attempts in determining stationarity: Initially I attempted to plot the distribution of amplitudes and phases to determine stationarity but I ran into problems here because I only have 50 samples per test case and half of my samples are split into two groups (read the unknown section). I tried to solve this issue by performing an RMS moving function to possibly normally distribute the values (This worked for taking a moving rms calculation with a length of 1 period (of the signal) and shifted by 1 sample). This somewhat worked for the distribution of phase but the variation and accuracy of my Anderson-Darling test would vary drastically on changes to the input values of the RMS moving function. I stop this RMS method in trying to resample my data + still had the split of distributed phase values, and I suspect this would add statistical bias.

After this I came to the conclusion that I needed more samples for each trial ( more than 50 ) for stationarity to be declared. So recently I have been trying to recreate a signal from the complex valued number of each FFT performed in 1 trial. I would add noise that is Gaussian but at what level of noise I am unsure. I would then perform a similar moving RMS function to the signal to sample into a normal distribution. Lastly I would perform t-test and F test for variance. I am unsure if this would be a correct way to address my overarching goal "when should we stop collecting data?".

Other stuff I tried:

• Box Muller Transform for the complex values. (ran into the problem of half my samples being distributed, what i can assume normally in to different regimes)
• Read a bit of "Time Series: Data Analysis and Theory" by David R. Brillinger (chapter 4: Stochastic Properties of Finite Fourier Transforms). This discussed how my fourier coefficients are normally distributed and would have a "Wishart Distribution". But again I suspect I wouldn’t have enough data to represent this per trial run and ... again 25 of my samples per trial are distributed by some amplitude/phase away.
• Attempted to make some inference about spectral density. Since my complex value of my FFT are known, I know the spectral density "value". Additionally, if the spectral density is known then for a stationary process the the auto-correlation function could be obtained. But I think I lack understanding here to even approach this one; I read a paper that discussed possible stochastic tests for stationarity that would involve some evolutionary spectral density function.
• Root tests for stationarity (Adjusted Dickey-Fuller test)

Additional Questions I have:

• Is what I’m trying to achieve pointless? I.e. I would be able to do a statistical test but the sample space I’m working with is too small to make any significant statistical decision? Basically asking if I should just throw in the towel with this.
• If it is correct to resample my data for what I described above (moving rms) into a normal distribution is it appropriate to perform a t-test for the sample mean of the most recent trial run and prior trial (10 min ago) to declare stationarity? I do have a similar question but for the F-Test and variance of the recent and prior trial.
• In relation to above, I suspect that only looking at the most recent and prior trial runs for stationarity may not be correct. Is there correct approach in dealing with "delayed" stochastic processes? Moreover, since I am not collecting data continuously would I need to use a different approach?

I appreciate you taking the time to read this. If anything a response on just say "hey you should read this, that person tried to do something similar" would so helpful. Don’t be afraid to send dense material. I have some publications and an extensive math/science background. I’m somewhat unaware with Python, MATLAB is my strong suit.

Answer 1

The "Allan Deviation" (or ADEV) is an excellent tool for determining how long we can consider the signal to be stationary. Existing free tools to facilitate its computation are Stable32 for non-Python users, and Allantools in Python.

ADEV is typically used to determine frequency stability, but can be equally used to determine interesting characteristics of any noise process. Given the OP is interested in the stationarity of a specific frequency, then as long as this is sampled with a source of higher stability (this is important!) then the duration for stationarity can be determined from the ADEV measurement directly.

Below is a plot of ADEV from www.everythingrf.com for the frequency out of an oscillator. The horizontal axis shows the averaging time $\tau$ and the vertical axis shows the resulting deviation which is a measure of the rms noise over that averaging interval (this is a simplistic explanation to provide a high level, I provide links below to other posts I have made here that give the more specific mathematical details).

We see up to 100 seconds of averaging time that the deviation is going down at a rate of $1/\sqrt{\tau}$, which is consistent with averaging of stationary white noise! In fact, in this part of the curve, the noise process is indeed "White FM" or white frequency modulation: The waveform as viewed as frequency fluctuations vs time appears as white noise. Those familiar with averaging over $N$ samples of stationary white noise may recognize the similarity here in that the standard deviation of such an averaging process goes down as $1/\sqrt{N}$. Shortly after 100 seconds this bottoms out which is the indication that if we were to observe (and average) the output frequency for any longer we would not get any further improvement in reducing the noise (this is exactly what occurs when a system is no longer stationary!). What we also see as we continue to average out beyond 250 seconds, that the noise in our result for longer averages actually gets worst.

For the OP to do this measurement properly, a more stable sampling clock will be needed, as that will act as the "reference" of what good is. A straightforward approach would be to determine the phase from the captured samples (for a single tone, instead of a full FFT, consider multiplying by a quadrature tone ($e^{-j\omega_c t}$ and then passing the result through a low pass filtering to get the in-phase (I) and quadrature (Q) components (this is basically a single bin of the DFT). For the I and Q signal the instantaneous phase for each sample can be determined (using ATAN2) which can then be input into either of the Allan Deviation tools I linked using 'phase' mode.

For additional reading including additional points on how this is used for the purpose the OP seeks, please see these additional posts:

How to interpret ADEV plot?

This one provides more details in the math behind the ADEV plot:

How to interpret Allan Deviation plot for gyroscope?

This one mentions using ADEV for a single DFT bin spefifically:

What determines the accuracy of the phase result in a DFT bin?