synthesize_formant_speech.py

speech.py

"""
Code converted by Keunwoo from https://ccrma.stanford.edu/~jos/fp/Formant_Filtering_Example.html (JOS)
"""

import numpy as np
import scipy.signal as signal
import librosa

def gen_speech(F, sig=None, filename='speech'):
    nsecs = len(F)
    R = np.exp(-np.pi * BW / fs)# pope radii
    theta = 2 * np.pi * F / fs # pole angles
    poles = R * np.exp(1j * theta)

    # B = 1
    A = np.real(np.poly(np.concatenate([poles, poles.conj()], axis=0)))
    B = np.zeros(A.shape)
    B[0] = 1
    r, p, f = signal.residuez(B, A)
    As = np.zeros((nsecs, 3), dtype=np.complex)
    Bs = np.zeros((nsecs, 3), dtype=np.complex)

    for idx, i in enumerate(range(1, 2*nsecs+1, 2)):
        j = i - 1

        Bs[idx] = [r[j] + r[j+1], -(r[j] * p[j+1] + r[j+1] * p[j]), 0]
        As[idx] = [1, -(p[j] + p[j+1]), p[j] * p[j+1]]

    sos = np.concatenate([As, Bs], axis=1)
    iperr = np.abs(np.imag(sos)) / (np.abs(sos) + 1e-10)
    sos = np.real(sos)
    Bh, Ah = signal.sos2tf(sos)
    nfft = 512

    H = np.zeros((nsec+1, nfft))
    for i in range(nsecs):
        Hiw, w = signal.freqz(Bs[i, :], As[i, :])
        H[i+1, :] = np.conj(Hiw[:])

    H[0, :] = np.sum(H[1:, :], axis=0)

    # gen signal if needed
    if sig is None:
        nsamps = 2048
        f0 = 200
        w0T = 2 * np.pi * f0 / fs
        nharm = np.floor((fs/2)/f0)

        sig = np.zeros(nsamps)
        nums = np.arange(0, nsamps)
        for i in np.arange(nharm):
            sig += np.cos(i * w0T * nums)

        sig = sig / np.max(np.abs(sig))
        
    speech = signal.lfilter([1], A, sig)
    speech = speech / np.max(np.abs(speech))
    librosa.output.write_wav(filename + '.wav', y=speech, sr=fs, norm=True)