Compute spectrum of sample-and-hold noise

SH1.py

# see https://lists.columbia.edu/pipermail/music-dsp/2015-November/000424.html
# psd derivation due to Ethan Duni
import numpy as np
from numpy.fft import fft, fftfreq
import matplotlib.pyplot as plt

N = 16384*2*2*2*2*2 # FFT size

y = (np.random.random(N) * 2) - 1 # ~U(-1,1)
r = np.random.random(N) # ~U(0,1)
x = np.empty(N) # (r[n]<P)?x[n-1]:y[n]

# hold process with probability of P for holding
P = 0.9
x[0] = 0
for i in range(1,N-1):
    if r[i] < P:
        x[i] = x[i-1]
    else:
        x[i] = y[i]
    
# FFT of x
spectrum = fft(x,N)
spectrum *= (1.0 / N)
magSqSpect = (np.abs(spectrum)) ** 2

# Verify that power of x and spectrum match:
meanTDPower = sum(x*x) / N
meanFDPower = sum(magSqSpect)
print meanFDPower / meanTDPower # 1.0

# Calculate psd. psd is in units of power / radian
w = fftfreq(N)*np.pi*2 # (w from 0 to 0.5 then -0.5 to 1
psd = ((1.0 - P**2)/(1.0 - 2*P*np.cos(w) + P**2))

psd = psd * (1.0/3.0) # scale normalized psd by power of y[n]

# magSqSpect was in units of power / bin-width
# bin width is 2pi/N radians.
# multiply by N/2pi to get power / radian to match psd units
magSqSpect = magSqSpect * N/(2*np.pi)

plt.plot(w, magSqSpect)
plt.plot(w, psd, color='red')
plt.xscale('log')
plt.yscale('log')
plt.show()