python语音信号处理_python实现语音信号处理常用度量方法
信噪比(SNR)
有用信号功率与噪声功率的比(此处功率为平均功率)，也等于幅度比的平方
其中：P_{signal}为信号功率(平均功率或者实际功率)；P_{noise}为噪声功率；A_{signal}为信号幅度；A_{noise}为噪声幅度值，功率等于幅度值的平方
MATLAB版本代码
#信号与噪声长度应该一样
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)
tensorflow版本SNR
deftf_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
defVolodymyr_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
Volodymyr Kuleshov论文实现方法
批注：这里的1e-6和1e-8，目的是为了防止出现Nan值，如果没有这个需求可以去除
numpy版本代码
defnumpy_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
峰值信噪比(PSNR)
表示信号的最大瞬时功率和噪声功率的比值，最大瞬时功率为语音数据中最大值得平方。
defpsnr(label, logits):
MAX= np.max(label) ** 2 #信号的最大平时功率
MSE = np.mean((label - logits) ** 2)return np.log10(MAX / MSE)
分段信噪比(SegSNR)
由于语音信号是一种缓慢变化的短时平稳信号，因而在不同时间段上的信噪比也应不一样。为了改善上面的问题，可以采用分段信噪比。分段信噪比即是先对语音进行分帧，然后对每一帧语音求信噪比，最好求均值。
MATLAB版本的代码
function [segSNR] =Evaluation(clean_speech,enhanced)
N= 25*16000/1000; %length of the segment interms 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-1sum1=0;
sum2=0;for n = mN +1 : mN+N
sum1= sum1 +clean_speech(n)^2;
sum2= sum2 +(enhanced{i}(n) - clean_speech(n))^2;
end
r= 10*log10(sum1/sum2);if r>55r= 55;
elseif r< -10r= -10;
end
segSNR(i)= segSNR(i) +r;
end
segSNR(i)= segSNR(i)/M;
end
View Code
python代码
defSegSNR(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 inrange(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)
对数拟然对比度(log Likelihood Ratio Measure)
坂仓距离测度是通过语音信号的线性预测分析来实现的。ISD基于两组线性预测参数(分别从原纯净语音和处理过的语音的同步帧得到)之间的差异。LLR可以看成一种坂仓距离(Itakura Distance,IS)但是IS距离需要考虑模型增益。而LLR不需要考虑模型争议引起的幅度位移，更重视整体谱包络的相似度。
PESQ
ITU-T的全系列参考目标语音质量测量系列
1997年的P.861(PSQM)
2001年的P.862(PESQ)，后来补充了P.862.1，P.862.2(宽带测量)，P.862.3(应用指南)
2011年的P.863(POLQA)
PESQ是用于语音质量评估的一种方法，ITU提供了C语言代码，下载请点击这里，但是在使用之前我们需要先编译C脚本，生成可执行文件exe
编译方式为：在命令行进入下载好的文件
cd \Software\source
gcc -o PESQ *.c
经过编译，会在当前文件夹生成一个pesq.exe的可执行文件
使用方式为：
命令行进入pesq.exe所在的文件夹
执行命令：pesq 采样率 "原始文件路径名" "劣化文件路径名”
回车
等待结果即可，值越大，质量越好。
例如：pesq +16000 raw.wav processed.wav
感知客观语音质量评估(POLQA)
POLQA是PESQ的继承者(ITU-T P.862建议书)。POLQA避免了当前P.862型号的弱点，并且扩展到处理更高带宽的音频信号。进一步的改进针对具有许多延迟变化的称为信号和信号的时间的处理。与P.862类似，POLQA支持普通电话频段(300-3400 Hz)的测量，但此外它还具有第二种操作模式，用于评估宽带和超宽带语音信号中的HD-Voice(50-14000)赫兹)。
POLQA是全参考算法，并且在对应的参考和测试信号的摘录的时间对准之后逐个样本地分析语音信号。POLQA可用于为网络提供端到端(E2E)质量评估，或表征各个网络组件。
POLQA结果主要是模型平均意见得分(MOS)，涵盖从1(差)到5(优秀)的范围。
对数谱距离(LSD)
对数谱距离Log Spectral Distance，LSD是两个频谱之间的距离度量。也称为“对数谱失真”
式中，l和m分别为频率索引和帧索引，M为语音帧数，L为频点数，\hat{S}(l, m)和S(l, m)分别为估计音频和宽带音频经过短时短时傅里叶变换后的频谱。
numpy版本
#方法一
defnumpy_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)))returnlsd#方法二
defget_power(x):
S= librosa.stft(x, n_fft=2048) #(1 + n_fft/2, n_frames)
S = np.log10(np.abs(S) ** 2)returnSdefcompute_log_distortion(labels, logits):"""labels和logits数据维度为 (batch_size, seq_len, 1)"""avg_lsd=0
batch_size=labels.shape[0]for i inrange(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+=lsdreturn avg_lsd / batch_size
tensorflow版本
defget_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.)
returnSdeftf_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版本的值和tensorflow版本的值一样，可以使用下面的代码
#numpy版本一：处理一个batch，(batch, seq_len, 1)
defnumpy_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) ** 2lsd= np.mean(np.sqrt(np.mean(original_target_squared, axis=0)))returnlsd#numpy版本二：处理一个batch，(batch, seq_len, 1)
defget_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)returnSdefcompute_log_distortion(labels, logits):
avg_lsd=0
batch_size=labels.shape[0]for i inrange(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+=lsdreturn avg_lsd /batch_size#tensorflow版本
defget_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.)returnSdeftf_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
View Code
批注：librosa.stft中center设为False，和np.log10中加1e-8，目的是为了最终的值和tensorflow版本的lsd值相近，如果没有这个需求可以去除。这里tf.log中加9.677e-9是为了和numpy中的值相近，如果没有这个需求可以去除
短时客观可懂度(STOI)
下载一个 pystoi 库：pip install pystoi
STOI 反映人类的听觉感知系统对语音可懂度的客观评价，STOI 值介于0~1 之间，值越大代表语音可懂度越高，越清晰。
from pystoi importstoi
stoi_score= stoi(label, logits, fs_sig=16000)
加权谱倾斜测度(WSS)
WSS值越小说明扭曲越少，越小越好，范围
参考文献
度量方法仓库