Illustration of Wiener-Khintchine

khintchine.py

import numpy as np
import matplotlib.pyplot as plt

# Signal frequencies, in Hz
f = np.array([10, 20, 40, 80])

# Sampling duration
T = 1.0

# Sampling frequency (in Hz)
fs = 200

# Number of samples
N = int(fs / T)

t = np.linspace(0, 1, N)

# Select which frequency to use at each time-step
s = (np.random.random(N) * len(f)).astype(int)

x = np.sin(f[s] * 2 * np.pi * t)
X = np.fft.fft(x)

x = x - np.mean(x)
XX = np.fft.ifft(np.correlate(x, x, mode='full'))

fig, (ax0, ax1, ax2) = plt.subplots(1, 3)
ax0.plot(t, x)

ax1.plot(np.linspace(-fs, fs, len(X)), np.fft.fftshift(X))
ax1.vlines(f, np.min(X), np.max(X), 'r', linewidth=3)

ax2.plot(np.linspace(-fs, fs, len(XX)), np.fft.fftshift(XX))
ax2.vlines(f, np.min(XX), np.max(XX), 'r', linewidth=3)

plt.show()