信号分析 / Signal Analysis

def spectrum(x, samprate):
    """
    compute amplitude spectrum
    Usage:
        freq, ampspec = spectrum(sig, samprate)
    """
    y = scipy.fft(x, nextpow2(len(x)))
    #y = scipy.fft(x)
    y_tmp = np.abs(y)
    y_ampspec = y_tmp[:len(y)//2+1]
    y_ampspec[1:(len(y)//2)] *= 2.0
    freq = scipy.linspace(start=0, stop=samprate/2, num=len(y)/2+1)
    return(freq, y_ampspec / len(x))

def plot_spectrum(x, samprate, xlim=None, ylim=None):
    f, ampspec = spectrum(x, samprate)
    plt.plot(f, ampspec)
    plt.xlabel("Frequency (Hz)")
    plt.ylabel("Amplitude Spectrum")
    plt.xlim(xlim)
    plt.ylim(ylim)
    plt.grid()
    plt.show()

def plot_magnitude(x, samprate, xlim=(20, 20000), ylim=(-100, 0)):
    f, spec = spectrum(x, samprate)
    mag = 10 * np.log10(spec ** 2)
    mag = mag - max(mag)
    plt.semilogx(f, mag)
    plt.xlabel("Frequency (Hz)")
    plt.ylabel("Power Spectrum (dB)")
    plt.xlim(xlim)
    plt.ylim(ylim)
    plt.grid()
    plt.show()

def spectral_centroid(x, samprate):
    freq, spec = spectrum(x, samprate)
    # remove DC and nyquist
    powspec = spec[1 : -1] ** 2
    f = freq[1 : -1]
    # compute spectral centroid (of log-frequency and power spectrum)
    if sum(powspec) != 0.0:
        speccent = 2**(sum(np.log2(f) * powspec) / sum(powspec))
    else:
        speccent = None
    return(speccent)

def spectral_centroid_linlin(x, samprate):
    freq, spec = spectrum(x, samprate)
    # remove DC and nyquist
    ampspec = spec[1 : -1]
    f = freq[1 : -1]
    # compute spectral centroid (of linear-frequency and amplitude spectrum)
    if sum(ampspec) != 0.0:
        speccent = sum(f * ampspec) / sum(ampspec)
    else:
        speccent = None
    return(speccent)

def spectral_centroid_moving(x, samprate, window_size=8192, hop_size=1024):
    """
    Calculate spectral centroid of a signal as in moving average.
    Usage:
        x, samprate = soundfile.read("audio.wav")
        sc, t = spectral_centroid_moving(x, samprate,
                                         window_size=8192, hop_size=1024)
        matplotlib.pyplot.plot(t/samprate, sc)
    """
    num_frames = (len(x) + window_size-1) // hop_size
    xx = np.concatenate((x, np.zeros(2 * window_size)))
    spec_cents = np.zeros(num_frames)
    for n in range(num_frames):
        xxx = xx[hop_size*n : hop_size*n+window_size-1]
        spec_cents[n] = spectral_centroid(xxx, samprate)
    t = np.arange(num_frames) * hop_size
    return(spec_cents, t)

def hilbert(x):
    """
    Hilbert Transform
    """
    # add zero at the end if signal length is odd
    if (len(x) % 2 == 1):
        flag = True
        x = np.concatenate((x, 0))
    else:
        flag = False
    # compute analytical signal
    xx = scipy.fft(x)
    hh = [1] + [2]*(len(xx)//2-1) + [1] + [0]*(len(xx)//2-1)
    yy = xx * hh
    y = scipy.ifft(yy).real
    # trim trailing zero to match the original signal length
    if (flag == True):
        return(y[:-1])
    else:
        return(y)

def envelope(x):
    """
    Signal Envelope using Hilbert Transform
    """
    return(np.abs(hilbert(x)))

def reverb_time_sabine(width, depth, height, alpha):
    """
    estimate reverberation time using Sabine's equation
    """
    S = (width * height + width * depth + depth * height) * 2.0
    V = width * depth * height
    return((0.161 * V) / (alpha * S))

def reverb_envelope(ir):
    """
    Schroeder's reverberation envelope
    """
    return(np.flipud(np.cumsum(np.flipud(ir ** 2))))

def reverb_time(ir, samprate, T=[-5,-35,-60]):
    """
    estimate reverberation time from impulse response
    """
    # compute reverberation envelope in dB
    t = np.linspace(0, len(ir) / samprate, len(ir))
    revenv = reverb_envelope(ir)
    irdb = 10 * scipy.log10(revenv)
    irdb = irdb - max(irdb)
    # linear fit between -5dB and -35dB and find the time of being -60dB
    ind = [sum(irdb >= T[0]), sum(irdb >= T[1])]
    p = np.polyfit(irdb[ind[0]:ind[1]], t[ind[0]:ind[1]], 1)
    rt = T[2] * p[0]
    # return reverb time
    return(rt)

def ambisonics_AtoB(A, B, C, D):
    """
    ambisonics A to B format conversion
    """
    W = ( A + B + C + D) / 2
    X = (-A + B + C - D) / 2   # front-back
    Y = ( A + B - C - D) / 2   # left-right
    Z = (-A + B - C + D) / 2   # up-down
    return(W, X, Y, Z)

##### SIGNAL PROCESSING #######################################################

def normalize(x):
    """
    Signal normalize to [-1, +1]
    """
    return(x / np.max(np.abs(x)))

def fadein(x, ms, samprate):
    """
    apply raised-cosine fade-in at the beginning of a signal
    """
    winsize = round(ms / 1000 * samprate)
    y = x
    if winsize != 0:
        t = np.arange(winsize) / winsize
        w = (np.cos((t+1) * math.pi) + 1) / 2
        y[:winsize] *= w
    return(y)

def fadeout(x, ms, samprate):
    """
    apply raised-cosine fade-out at the end of a signal
    """
    winsize = round(ms / 1000 * samprate)
    y = x
    if winsize != 0:
        t = np.arange(winsize) / winsize
        w = (np.cos(t * math.pi) + 1) / 2
        y[-winsize:] *= w
    return(y)

def freeze(x):
    """
    Freezing audio by randomizing phase component.
    Usage:
        filename = 'in.wav'
        x, samprate = soundfile.read(filename)
        sinfo = soundfile.info(filename)
        x2 = freeze(x)
        soundfile.write("out.wav", x2, sinfo.samplerate, sinfo.subtype)
    """
    # FFT
    y = scipy.fft(x)
    mag = scipy.absolute(y)
    phase = scipy.angle(y)
    # phase randomization
    p2 = phase[1:len(phase)//2]
    np.random.seed(0)
    perm = np.random.permutation(len(p2))
    p2p = p2[perm]
    phaseR = scipy.r_[phase[0], p2p, phase[len(phase)//2], -p2p[::-1]]
    # IFFT
    y2 = mag * scipy.exp(1j * phaseR)
    x2 = scipy.real(scipy.ifft(y2))
    return(x2)