I am referring to the work of Stephen A. Billings on "Identification of a class of nonlinear systems using correlation analysis" from the year 1978, where it is mentioned that the $i^{\text{th}}$ dimensional autocorrelation function for a zero mean white Gaussian process $x(t)$ with spectral density of 1 W/cycle is as follows: \begin{equation} \overline{x(t_1)x(t_2)\dots x(t_i)} = \begin{cases}0,& \text{odd } i\\ \sum_i{ \prod_{n\ne m} \delta(t_n-t_m)},& \text{even } i\\ \end{cases} \end{equation}

I would like to use this result in my work and check if this is true. How do I create a signal with prescribed spectral density in MATLAB? and how do I perform autocorrelation for $i > 2$ in MATLAB?


1 Answer 1


The statement is provably true:

  1. $x(t)$ is zero mean and Gaussian with variance 1.
  2. $$E[x(t_i)x(t_j)] = \begin{cases} 0 \mbox{ for } t_i \not= t_j \\ 1 \mbox{ for } t_i = t_j\end{cases}$$

For $$E[ x(t_1) x(t_2) \ldots x(t_N)]$$

  1. For $N$ odd: There will always be a "left over" $E[x(t_N)] = 0$ multiplying the result, so it's always zero.

  2. For $N$ even: There will be a non-zero value only if $t_i = t_j$ for distinct $i$ and $j$. For the specific case where all times are equal,


where !! represents the double factorial function (not two applications of the factorial). For each pair of $t_i$, $t_j$ that are equal, it will generate another 1 (since $\sigma = 1$ here).

Note that the paper uses time averages and the above uses ensemble averages. The paper assumes ergodicity of the time series, so the two are equivalent.

Much as it goes against my better judgement, here's some code that does this. The output validates the "proof":

For 2: should be 1 Calculated: 0.9946066112024137 
For 3: should be 0 Calculated: 0.014338673370010435 
For 4: should be 3 Calculated: 3.0088602833619134 
For 5: should be 0 Calculated: 0.12289469288233784 
For 6: should be 15 Calculated: 15.20500109670608 
For 7: should be 0 Calculated: 1.6572818702690346

Code Below

import numpy as np
import math
from matplotlib import pyplot as plt

N = 1000
Nruns = 100000

def do_simulation(M):
    indices = 11*np.ones(M, dtype=int)
    sum_product = 0
    for idx in np.arange(Nruns):
        x = np.random.normal(size=N)
        sum_product = sum_product + math.prod(x[indices])

    return sum_product/Nruns

def number_should_be(M):
    if (M % 2 == 0):
        return math.prod(range(M-1, 0, -2))
        return 0

print("For 2: should be " + str(number_should_be(2)) + " Calculated: " + str(do_simulation(2)))
print("For 3: should be " + str(number_should_be(3)) + " Calculated: " + str(do_simulation(3)))
print("For 4: should be " + str(number_should_be(4)) + " Calculated: " + str(do_simulation(4)))
print("For 5: should be " + str(number_should_be(5)) + " Calculated: " + str(do_simulation(5)))
print("For 6: should be " + str(number_should_be(6)) + " Calculated: " + str(do_simulation(6)))
print("For 7: should be " + str(number_should_be(7)) + " Calculated: " + str(do_simulation(7)))
  • $\begingroup$ Thank you for the detailed explanation. Do you have an answer why the spectral density of the signal should be 1 W/cycle? And How is this achieved in your Python code? Does the random.normal function automatically returns a signal with spectral density 1? $\endgroup$
    – Neuling
    Commented Jan 16 at 12:17
  • $\begingroup$ Also, I would like to ask whether you are acquainted with this paper? The paper is pretty old, some of the references are not available, and the explanation in the paper is not quite comprehensible. I have a Wiener-Hammerstein system (dynamic linear block followed by static nonlinear block and by dynamic linear block) with input $u$ and output $y$. I would like to understand, how the parameters of these individual blocks are estimated in this paper. $\endgroup$
    – Neuling
    Commented Jan 16 at 14:42
  • $\begingroup$ @Neuling I'm just assuming that the 1 W/cycle means a standard deviation of 1. That's achieved by np.random.normal(size=N) which has a mean of 0 and a standard deviation of 1 as the default values. Your second comment: that sounds like a whole other question. $\endgroup$
    – Peter K.
    Commented Jan 21 at 15:43

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.