python實現語音訊號處理常用度量方法
- 2019 年 10 月 3 日
- 筆記
信噪比(SNR)
有用訊號功率與雜訊功率的比（此處功率為平均功率），也等於幅度比的平方
$$SNR(dB)=10log_{10}frac{sum_{n=0}^{N-1}s^2(n)}{sum_{n=0}^{N-1}d^2(n)}=10*log_{10}(frac{P_{signal}}{P_{noise}})=20*log_{10}(frac{A_{signal}}{A_{noise}})$$
$$SNR(dB)=10log_{10}frac{sum_{n=0}^{N-1}s^2(n)}{sum_{n=0}^{N-1}[x(n)-s(n)^2]}$$
其中：
$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)
View Code
python程式碼
def numpy_SNR(origianl_waveform, target_waveform): # 單位 dB signal = np.sum(origianl_waveform ** 2) noise = np.sum((origianl_waveform - target_waveform) ** 2) snr = 10 * np.log10(signal / noise) return snr
$$np.linalg.norm(x)=sqrt{x_1^2+x_2^2+…+x_n^2}$$
這個公式和上面是一樣的
def wav_snr(ref_wav, in_wav):# 如果ref wav稍長，則用0填充in_wav if (abs(in_wav.shape[0] - ref_wav.shape[0]) < 10): pad_width = ref_wav.shape[0] - in_wav.shape[0] in_wav = np.pad(in_wav, (0, pad_width), 'constant') else: print("錯誤：參考wav與輸入wav的長度明顯不同") return -1 # 計算 SNR norm_diff = np.square(np.linalg.norm(in_wav - ref_wav)) if (norm_diff == 0): print("錯誤：參考wav與輸入wav相同") return -1 ref_norm = np.square(np.linalg.norm(ref_wav)) snr = 10 * np.log10(ref_norm / norm_diff) return snr
峰值信噪比(PSNR)
表示訊號的最大瞬時功率和雜訊功率的比值，最大瞬時功率為語音數據中最大值得平方。
$$SNR(dB)=10*log _{10}(frac{MAX(P_{signal})}{P_{noise}})=10log_{10}frac{MAX[s(n)]^2}{d^2(n)}$$
$$SNR(dB)=10log_{10}frac{MAX[s(n)]^2}{frac{1}{N}sum_{n=0}^{N-1}[x(n)-s(n)]^2}=20log_{10}frac{MAX[s(n)]}{sqrt{MSE}}$$
import numpy as np def psnr(ref_wav, in_wav): MSE = numpy.mean((ref_wav - in_wav) ** 2) MAX = np.max(ref_wav) # 訊號的最大平時功率 return 20 * np.log10(MAX / np.sqrt(MSE))
分段信噪比(SegSNR)
由於語音訊號是一種緩慢變化的短時平穩訊號，因而在不同時間段上的信噪比也應不一樣。為了改善上面的問題，可以採用分段信噪比。分段信噪比即是先對語音進行分幀，然後對每一幀語音求信噪比，最好求均值。
MATLAB版本的程式碼
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
View Code
python程式碼
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 # num_frame = (len(ref_wav)-windowsize+shift) // shift num_frame = (len(ref_wav) - windowsize) // shift + 1 # 計算幀的數量 SegSNR = np.zeros(num_frame) # 計算每一幀的信噪比 for i in range(0, 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
PESQ是用於語音品質評估的一種方法，ITU提供了C語言程式碼，下載請點擊這裡，但是在使用之前我們需要先編譯C腳本，生成可執行文件exe
編譯方式為：在命令行進入下載好的文件
- cd Softwaresource
- gcc -o PESQ *.c
經過編譯，會在當前文件夾生成一個pesq.exe的可執行文件
使用方式為：
- 命令行進入pesq.exe所在的文件夾
- 執行命令：pesq 取樣率 “原始文件路徑名” “劣化文件路徑名”
- 回車
- 等待結果即可，值越大，品質越好。
- 例如：pesq +16000 raw.wav processed.wav
對數譜距離(Log Spectral Distance)
對數譜距離Log Spectral Distance是兩個頻譜之間的距離度量（用分貝表示）。兩個頻譜$P(W)$和$hat{P}(w)$之間的對數譜距離被定義為：
$$D_{LS}=sqrt{frac{1}{2pi}int_{-pi}^{pi}[10*log _{10}frac{P(w)}{hat{P}(w)}]^2dw}$$
其中，$p(w)$和$hat{P}(w)$是功率譜。對數譜距離是時多對稱的。
def numpy_LSD(origianl_waveform, target_waveform): """ 比較原始和目標音頻之間的對數譜距離（LSD），也稱為對數譜失真， 是兩個頻譜之間的距離測量值（以dB表示） """ print("數據形狀為", origianl_waveform.shape) print("數據類型為", type(origianl_waveform)) original_spectrogram = librosa.core.stft(origianl_waveform, n_fft=2048) target_spectrogram = librosa.core.stft(target_waveform, n_fft=2048) original_log = np.log10(np.abs(original_spectrogram) ** 2) target_log = np.log10(np.abs(target_spectrogram) ** 2) original_target_squared = (original_log - target_log) ** 2 target_lsd = np.mean(np.sqrt(np.mean(original_target_squared, axis=0))) return target_lsd
參考文獻：
