def histo(fn):
    import scikits.audiolab as au

    (sp, sf, b) = au.wavread(fn)

    sp = sp.conj().transpose()
    sfmmean, nmeas = 0, 0
    his=zeros((65536,), dtype=numpy.int)
    windd = hanning(2048).conj().transpose()

    for i in range(len(sp[1])):
        for j in (0, 1):
            t = round(sp[j, i] * 32768 + 32769)
            his[t] += 1


            if i % 1024 == 0 and i > 1024:
                w = sp[0, i-2048:i]
                w = w * windd
                wt = pylab.fft(w)
                wt = wt * conj(wt)
                wt = maximum(wt, .0000000001)
                am = sum(wt) / 2048
                gm = exp(sum(log(wt)) / 2048)
                sfm = gm / am
                if am > .5:
                    nmeas += 1
                    sfm = 10 * log10(sfm)
                    sfmmean += sfm

    his = maximum(his/float(sum(his)), .000000000001)
    xax = arange(-32768,32768)
    semilogy(xax,his)