Skip to content

Instantly share code, notes, and snippets.

@RossBencina
Last active November 5, 2015 14:52
Show Gist options
  • Save RossBencina/a15a696adf0232c73a55 to your computer and use it in GitHub Desktop.
Save RossBencina/a15a696adf0232c73a55 to your computer and use it in GitHub Desktop.
Compute spectrum of sample-and-hold noise
# 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()
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment