python实现语音常用度量方法的代码详解

# 信号与噪声长度应该一样
function snr=SNR_singlech(Signal,Noise)

P_signal = sum(Signal-mean(Signal)).^2;     # 信号的能量
P_noise = sum(Noise-mean(Noise)).^2;     # 噪声的能量
snr = 10 * log10(P_signal/P_noise)

def tf_compute_snr(labels, logits):
    # labels和logits都是三维数组 (batch_size, wav_data, 1)
    signal = tf.reduce_mean(labels ** 2, axis=[1, 2])
    noise = tf.reduce_mean((logits - labels) ** 2, axis=[1, 2])
    noise = tf.reduce_mean((logits - labels) ** 2 + 1e-6, axis=[1, 2])
    snr = 10 * tf.log(signal / noise) / tf.log(10.)
    # snr = 10 * tf.log(signal / noise + 1e-8) / tf.log(10.)
    snr = tf.reduce_mean(snr, axis=0)
    return snr

def Volodymyr_snr(labels, logits):
    # labels和logits都是三维数组 (batch_size, wav_data, 1)
    noise = tf.sqrt(tf.reduce_mean((logits - labels) ** 2 + 1e-6, axis=[1, 2]))
    signal = tf.sqrt(tf.reduce_mean(labels ** 2, axis=[1, 2]))
    snr = 20 * tf.log(signal / noise + 1e-8) / tf.log(10.)
    avg_snr = tf.reduce_mean(snr, axis=0)
    return avg_snr

def numpy_SNR(labels, logits):
    # origianl_waveform和target_waveform都是一维数组 (seq_len, )
    # np.sum实际功率;np.mean平均功率，二者结果一样
    signal = np.sum(labels ** 2)
    noise = np.sum((labels - logits) ** 2)
    snr = 10 * np.log10(signal / noise)
    return snr

def psnr(label, logits):
    MAX = np.max(label) ** 2  # 信号的最大平时功率
    MSE = np.mean((label - logits) ** 2)
    return np.log10(MAX / MSE)

function [segSNR] = Evaluation(clean_speech,enhanced)

N = 25*16000/1000; %length of the segment in terms of samples
M = fix(size(clean_speech,1)/N); %number of segments
segSNR = zeros(size(enhanced));
for i = 1:size(enhanced,1)
    for m = 0:M-1
        sum1 =0;
        sum2 =0;
        for n = m*N +1 : m*N+N
            sum1 = sum1 +clean_speech(n)^2;
            sum2 = sum2 +(enhanced{i}(n) - clean_speech(n))^2;
        end
        r = 10*log10(sum1/sum2);
        if r>55
            r = 55;
        elseif r  -10
            r = -10;
        end
       
        segSNR(i) = segSNR(i) +r;
    end
    segSNR(i) = segSNR(i)/M;
end

def SegSNR(ref_wav, in_wav, windowsize, shift):
    if len(ref_wav) == len(in_wav):
        pass
    else:
        print('音频的长度不相等!')
        minlenth = min(len(ref_wav), len(in_wav))
        ref_wav = ref_wav[: minlenth]
        in_wav = in_wav[: minlenth]
    # 每帧语音中有重叠部分，除了重叠部分都是帧移，overlap=windowsize-shift
    # num_frame = (len(ref_wav)-overlap) // shift
    #           = (len(ref_wav)-windowsize+shift) // shift
    num_frame = (len(ref_wav) - windowsize + shift) // shift  # 计算帧的数量

    SegSNR = np.zeros(num_frame)
    # 计算每一帧的信噪比
    for i in range(num_frame):
        noise_frame_energy = np.sum(ref_wav[i * shift: i * shift + windowsize] ** 2)  # 每一帧噪声的功率
        speech_frame_energy = np.sum(in_wav[i * shift: i * shift + windowsize] ** 2)  # 每一帧信号的功率
        SegSNR[i] = np.log10(speech_frame_energy / noise_frame_energy)

    return 10 * np.mean(SegSNR)

def SER(near_speech, far_echo):
    """signal to echo ratio, 信号回声比
    :param near_speech: 近端语音
    :param far_echo: 远端回声 
    """
    return 10*np.log10(np.mean(near_speech**2)/np.mean(far_echo**2))

def compute_ERLE(mic_wav, predict_near_end_wav):
    """ 
    :param mic_wav: 麦克风信号(y) = 近端语音(s) + 远端语音回声(s) + 噪声(v)
    :param predict_near_end_wav: 估计的近端语音信号 \hat{s}
    麦克风信号
    """
    mic_mear = np.mean(mic_wav**2)
    predict_near_end_wav = np.mean(predict_near_end_wav**2)
    ERLE = 10 * np.log10(mic_mear/predict_near_end_wav)
    return ERLE

from scipy.io import wavfile
from pesq import pesq

rate, ref = wavfile.read("./audio/speech.wav")
rate, deg = wavfile.read("./audio/speech_bab_0dB.wav")

print(pesq(rate, ref, deg, 'wb'))
print(pesq(rate, ref, deg, 'nb'))

def pesq2mos(pesq):
    """ 将PESQ值[-0.5, 4.5]映射到MOS-LQO得分[1, 4.5]上，映射函数来源于：P.862.1 """
    return 0.999 + (4.999 - 0.999) / (1 + np.exp(-1.4945 * pesq + 4.6607))

def mos2pesq(mos):
    """ 将MOS-LQO得分[1, 4.5]映射到PESQ值[-0.5, 4.5]上，映射函数来源于：P.862.1"""
    inlog = (4.999 - mos) / (mos - 0.999)
    return (4.6607 - np.log(inlog)) / 1.4945

# 方法一
def numpy_LSD(labels, logits):
    """ labels 和 logits 是一维数据 (seq_len,)"""
    labels_spectrogram = librosa.stft(labels, n_fft=2048)  # (1 + n_fft/2, n_frames)
    logits_spectrogram = librosa.stft(logits, n_fft=2048)  # (1 + n_fft/2, n_frames)

    labels_log = np.log10(np.abs(labels_spectrogram) ** 2)
    logits_log = np.log10(np.abs(logits_spectrogram) ** 2)
    # 先处理频率维度
    lsd = np.mean(np.sqrt(np.mean((labels_log - logits_log) ** 2, axis=0)))

    return lsd


# 方法二
def get_power(x):
    S = librosa.stft(x, n_fft=2048)  # (1 + n_fft/2, n_frames)
    S = np.log10(np.abs(S) ** 2)
    return S


def compute_log_distortion(labels, logits):
    """labels和logits数据维度为 (batch_size, seq_len, 1)"""
    avg_lsd = 0
    batch_size = labels.shape[0]
    for i in range(batch_size):
        S1 = get_power(labels[i].flatten())
        S2 = get_power(logits[i].flatten())
        # 先处理频率轴，后处理时间轴
        lsd = np.mean(np.sqrt(np.mean((S1 - S2) ** 2, axis=0)), axis=0)
        avg_lsd += lsd
    return avg_lsd / batch_size

def get_power(x):
    x = tf.squeeze(x, axis=2)  # 去掉位置索引为2维数为1的维度 （batch_size, input_size）
    S = tf.signal.stft(x, frame_length=2048, frame_step=512, fft_length=2048,
                       window_fn=tf.signal.hann_window)
    # [..., frames, fft_unique_bins]
    S = tf.log(tf.abs(S) ** 2) / tf.log(10.)
    # S = tf.log(tf.abs(S) ** 2 + 9.677e-9) / tf.log(10.)
    return S


def tf_compute_log_distortion(labels, logits):
    """labels和logits都是三维数组 (batch_size, input_size, 1)"""
    S1 = get_power(labels)  # [..., frames, fft_unique_bins]
    S2 = get_power(logits)  # [..., frames, fft_unique_bins]
    # 先处理频率维度，后处理时间维度
    lsd = tf.reduce_mean(tf.sqrt(tf.reduce_mean((S1 - S2) ** 2, axis=2)), axis=1)
    lsd = tf.reduce_mean(lsd, axis=0)
    return lsd

# numpy版本一：处理一个batch，(batch, seq_len, 1)
def numpy_LSD(labels, logits):
    """ labels 和 logits 是一维数据"""
    labels_spectrogram = librosa.stft(labels, n_fft=2048, hop_length=512, win_length=2048,
                                      window="hann", center=False)  # (1 + n_fft/2, n_frames)
    logits_spectrogram = librosa.stft(logits, n_fft=2048, hop_length=512, win_length=2048,
                                      window="hann", center=False)  # (1 + n_fft/2, n_frames)

    labels_log = np.log10(np.abs(labels_spectrogram) ** 2 + 1e-8)
    logits_log = np.log10(np.abs(logits_spectrogram) ** 2 + 1e-8)
    original_target_squared = (labels_log - logits_log) ** 2
    lsd = np.mean(np.sqrt(np.mean(original_target_squared, axis=0)))

    return lsd


# numpy版本二：处理一个batch，(batch, seq_len, 1)
def get_power1(x):
    S = librosa.stft(x, n_fft=2048, hop_length=512, win_length=2048,
                     window="hann", center=False)  # (1 + n_fft/2, n_frames)
    S = np.log10(np.abs(S) ** 2 + 1e-8)
    return S


def compute_log_distortion(labels, logits):
    avg_lsd = 0
    batch_size = labels.shape[0]
    for i in range(batch_size):
        S1 = get_power1(labels[i].flatten())
        S2 = get_power1(logits[i].flatten())
        # 先处理频率轴，后处理时间轴
        lsd = np.mean(np.sqrt(np.mean((S1 - S2) ** 2, axis=0)), axis=0)
        avg_lsd += lsd
    return avg_lsd / batch_size


# tensorflow版本
def get_power(x):
    x = tf.squeeze(x, axis=2)  # 去掉位置索引为2维数为1的维度 （batch_size, input_size）
    S = tf.signal.stft(x, frame_length=2048, frame_step=512, fft_length=2048,
                       window_fn=tf.signal.hann_window)
    # [..., frames, fft_unique_bins]
    S = tf.log(tf.abs(S) ** 2 + 9.677e-9) / tf.log(10.)
    return S


def tf_compute_log_distortion(labels, logits):
    # labels和logits都是三维数组 (batch_size, input_size, 1)
    S1 = get_power(labels)  # [..., frames, fft_unique_bins]
    S2 = get_power(logits)  # [..., frames, fft_unique_bins]
    # 先处理频率维度，后处理时间维度
    lsd = tf.reduce_mean(tf.sqrt(tf.reduce_mean((S1 - S2) ** 2, axis=2)), axis=1)
    lsd = tf.reduce_mean(lsd, axis=0)
    return lsd

from pystoi import stoi

stoi_score = stoi(label, logits, fs_sig=16000)