def many_psds(k=2,fs=1.0, b0=1.0, N=1024):
    """ compute average of many PSDs """
    psd=[]
    for j in range(k):
        print j
        x = noise.white(N=2*4096,b0=b0,fs=fs)
        f, tmp = noise.numpy_psd(x,fs)
        if j==0:
            psd = tmp
        else:
            psd = psd + tmp
    return f, psd/k
compute average of many PSDs