def ms(x):
    return (np.abs(x)**2.0).mean()

def normalize(y, x=None):
    if x is not None:
        x = ms(x)
    else:
        x = 1.0
    return y * np.sqrt( x / ms(y) )

def brown(N, state=None):
    state = np.random.RandomState() if state is None else state
    uneven = N%2
    X = state.randn(N//2+1+uneven) + 1j * state.randn(N//2+1+uneven)
    S = (np.arange(len(X))+1)# Filter
    y = (np.fft.irfft(X/S)).real
    if uneven:
        y = y[:-1]
    return normalize(y)

inst_freq = np.cumsum(brown(100))
inst_freq -= inst_freq.min() - 1
inst_freq /= inst_freq.max()
rectified_inst_freq = 1.0 / (np.floor(1.0 / inst_freq))

plt.plot(inst_freq)
plt.plot(rectified_inst_freq, 'r')
plt.show()
plt.plot(np.cos(np.cumsum(np.pi * inst_freq)))
plt.show()
plt.plot(np.cos(np.cumsum(np.pi * rectified_inst_freq)))
plt.show()