语音信号处理基础知识
1.均值
1.1.均值公式定义
x ‾ = x 1 + x 2 + ⋯ + x n n = ∑ j = 1 n x j n \overline{x} = \frac{x_1 + x_2 + \cdots + x_n}{n} = \frac{\sum^{n}_{j=1}x_j}{n} x=nx1+x2+⋯+xn=n∑j=1nxj
1.2.物理意义
一段语音信号的均值代表了语音的直流分量, 正常情况下音频系统的输入输出不存在交流分量。 如果发生了直流偏执的现象, 需要在信号预处理阶段去除直流分量。
1.3.代码实现
读取一段新闻语音文件, 计算它的均值。
from scipy.io import wavfile
import numpy as np
framerate, data = wavfile.read('./source.wav')
data_len = len(data)
sum = 0
for i in range(data_len):
sum += data[i]
mean_data = sum / data_len
print('mean of audio : %f' % mean_data)
mean of audio : -0.000320
numpy提供average函数来计算均值, 效果一致。
np.average(data)
-0.000320226621917049
可以看出， 均值接近0.
2.方差
2.1.方差的公式定义
s 2 = 1 n ∑ i = 1 n ( x i − x ‾ ) 2 s^2 = \frac{1}{n}\sum_{i=1}^{n}(x_i - \overline{x})^2 s2=n1i=1∑n(xi−x)2
2.2.物理意义
方差描述语音信号的波动范围, 交流分量的强弱, 即交流信号的平均功率。
2.3.代码实现
sum = 0
for i in range(data_len):
sum += (data[i] - mean_data) ** 2
var_data = sum / data_len
print('var of audio : %f' % var_data)
var of audio : 3197990.844967
numpy提供var函数用以计算方差, 效果一致。
np.var(data)
3197990.8449671054
3.均方差
3.1.均方差的公式定义
σ = 1 n ∑ i = 1 n ( x i − x ‾ ) 2 σ = \sqrt{\frac{1}{n}\sum_{i=1}^{n}(x_i - \overline{x})^2} σ=n1i=1∑n(xi−x)2
3.2.物理意义
均方差等于方差的方根, 和方差一样可以反映语音信号的离散程度。由于方差的平方计算造成量纲的倍数变化无法直观反映出偏离程度, 均方差的意义更为直观。
3.3.代码实现
sum = 0
for i in range(data_len):
sum += (data[i] - mean_data) ** 2
rmse_data = np.sqrt(sum / data_len)
print('rmse of data : %f' % rmse_data)
rmse of data : 1788.292718
numpy提供std函数用以计算均方差, 效果一致。
np.std(data)
1788.2927179203928
4.协方差
4.1.公式定义
C o v ( X , Y ) = ∑ i = 0 n ( X − X ‾ ) ( Y − Y ‾ ) n − 1 Cov(X, Y) = \frac{\sum_{i=0}^{n}(X-\overline{X})(Y - \overline{Y})}{n-1} Cov(X,Y)=n−1∑i=0n(X−X)(Y−Y)
4.2.物理意义
描述两个变量之间的变化趋势相关性
4.3.代码仿真理解
计算初始相位0， π之间
def cov(x, y):
avr_x = np.average(x)
avr_y = np.average(y)
sum = 0
for i in range(len(x)):
sum += (x[i] - avr_x) * (y[i] - avr_y)
return sum / (len(x) - 1)
import matplotlib.pyplot as plt
x = np.arange(0, 10 * np.pi, 0.1)
y_1 = np.sin(x)
plt.plot(x, y_1)
plt.xlabel('sample (n)')
plt.ylabel('amp')
plt.show()
初始相位为0的正弦信号图像
y_2 = np.sin(x + np.pi)
plt.plot(x, y_2)
[<matplotlib.lines.Line2D at 0x7efef0fbd2b0>]
初始相位为0的正弦信号图像,可以看出相位相差π的两个序列的变化趋势完全相反。计算这两个信号序列的协方差
print(cov(y_1, y_2))
-0.5002531220184325
numpy提供cov函数来计算协方差矩阵, 矩阵的次对角线是两个序列的协方差。
np.cov(y_1, y_2)
array([[ 0.50025312, -0.50025312],
[-0.50025312, 0.50025312]])
计算不同初始相位下, 两个正弦序列的协方差变化情况。
phi_list = np.linspace(0, 2 * np.pi, 100)
cov_list = [cov(np.sin(x), np.sin(x + phi)) for phi in phi_list]
plt.plot(phi_list, cov_list)
[<matplotlib.lines.Line2D at 0x7efef0af0208>]
5.分帧加窗以及频谱泄露
5.1.频谱泄露
计算离散傅里叶变换的前提是假设输入有限长度的信号序列是一个周期信号的一个完整周期, 周期长度为该序列的长度。离散傅里叶变换本质是以该序列长度为周期进行周期延拓后的周期信号, 计算这个无限周期序列的频谱图。由于dft计算出的频谱对应的是周期延拓后的周期信号, 如果有限长的信号的首尾不相等, 会使得周期延拓后的周期信号不连贯, 会带来一些本来不属于信号本身的频率分量。
下面定义dft实现如下:
import math
def dft(x):
""" :bref 计算输入序列的离散傅里叶变换 :param x : 输入序列 :param N : dft长度N :return spectrum : 频谱 """
N = len(x)
spec = np.zeros(N)
for k in range(N):
sum = 0
for n in range(N):
sum += x[n] * np.exp(-1j * 2 * np.pi * n * k / N)
spec[k] = abs(sum / N)
return spec
定义计算频谱对应的频率值函数:
def dft_freq(N, sample_rate):
""" :bref 计算对应的频域分辨率 :param N : DFT点数 :param sample_rate : 采样率 :return freqs : 对应的频域的采样率 """
freqs = sample_rate / N * np.array(range(N))
return freqs
采用频率为10的正弦信号, 使用五个整数倍周期作为信号的长度; 从绘制的波形图可以看出整数倍周期的首尾是连贯的。
def create_sin(Ts, freq, times):
sample_rate = 1 / Ts
freq = 10
T = 1 / freq
duration = times * T
samples = math.ceil(duration / Ts)
t = Ts * np.array(range(samples))
x = np.cos(2 * np.pi * freq * t)
return t, x
Ts = 0.01
t, x = create_sin(Ts, 10, 5)
plt.plot(t, x)
[<matplotlib.lines.Line2D at 0x7efeef8ec5c0>]
将该信号进行周期延拓三次, 可以看出周期延拓后的信号和原本一致。
def sig_append(x, times, Ts):
""" :bref 周期延拓 :param x : 需要周期延拓的信号 :param times : 延拓次数 :param Ts : 采样周期 :return t_add : 延拓后序列的时间 :return x_add : 延拓后的信号幅值 """
N = x.shape[0]
t_add = Ts * np.array(range(N * times))
x_add = np.zeros(N * times)
for i in range(times):
x_add[i * N:(i + 1) * N] = x[:]
return t_add, x_add
t_add, x_add = sig_append(x, 3, Ts)
plt.plot(t_add, x_add)
[<matplotlib.lines.Line2D at 0x7efeef56d860>]
使用dft计算该信号频谱, 可以看出频谱分量只有10Hz。
freq_dft = dft_freq(x.shape[0], 1 / Ts)
plt.plot(freq_dft, dft(x))
plt.xlabel('Freq (Hz)')
plt.xlim([0, freq_dft[-1] / 2])
plt.show()
将信号长度改为非整数倍周期5.2, 来看看时域频域的效果。
t_5_2, x_5_2 = create_sin(0.01, 10, 5.2)
plt.plot(t_5_2, x_5_2)
[<matplotlib.lines.Line2D at 0x7efeec7f6908>]
同样的, 将周期进行三倍延拓看看效果。可以发现在延拓交界点存在不期望的波形。
t_add, x_add = sig_append(x_5_2, 3, Ts)
plt.plot(t_add, x_add)
[<matplotlib.lines.Line2D at 0x7efeed7d9cf8>]
来看看频域的效果, 可以发现在10Hz分量之外还有其他的频率分量, 整个信号的能量从基频泄露了一部分到其他频率上, 这就是频谱泄露。
freq_dft = dft_freq(x_10_1.shape[0], 1 / Ts)
plt.plot(freq_dft, dft(x_10_1))
plt.xlabel('Freq (Hz)')
plt.xlim([0, freq_dft[-1] / 2])
plt.show()
可以通过对信号的首尾进行抑制来减弱首尾不连贯带来的频谱泄露, 即就是加窗。
汉宁窗的数学定义是:
w ( n ) = 0.5 × ( 1 − c o s ( 2 π n / ( N − 1 ) ) ) , 0 ≤ n ≤ N − 1 w(n) = 0.5 × (1 - cos(2\pi n / (N - 1))), 0 \le n \le N - 1 w(n)=0.5×(1−cos(2πn/(N−1))),0≤n≤N−1
定义生成汉宁窗的函数如下:
def create_hanning(N):
window = np.arange(N)
window = (1 - np.cos(window / (N - 1) * 2 * np.pi) + 1) / 2
return window
对5.2倍周期的信号进行汉宁窗加权, 可以看到首尾的信号被衰减了。
window = create_hanning(x_5_2.shape[0])
x_5_2_window = x_5_2 * window
plt.figure(0)
plt.title('hanning window')
plt.plot(range(window.shape[0]), window)
plt.figure(1)
plt.title('after window')
plt.plot(t_5_2, x_5_2_window)
[<matplotlib.lines.Line2D at 0x7efeec381358>]
绘制对应的频谱图, 可以看出非10Hz的信号分量少了很多。
freq_dft = dft_freq(x_5_2_window.shape[0], 1 / Ts)
plt.plot(freq_dft, dft(x_5_2_window))
plt.xlabel('Freq (Hz)')
plt.xlim([0, freq_dft[-1] / 2])
plt.show()
def frame_sig(sig,frame_len,frame_step,winfunc=lambda x:numpy.ones((x,))):
""" :bref 将信号进行分帧加窗处理 :param sig: 需要分帧加窗的语音信号. :param frame_len: 帧长. :param frame_step: 帧移. :param winfunc: 窗函数, 默认为矩形窗. :returns: 分帧后的数据. """
slen = len(sig)
frame_len = int(round_half_up(frame_len))
frame_step = int(round_half_up(frame_step))
if slen <= frame_len:
numframes = 1
else:
numframes = 1 + int(math.ceil((1.0*slen - frame_len)/frame_step))
padlen = int((numframes-1)*frame_step + frame_len)
zeros = numpy.zeros((padlen - slen,))
padsignal = numpy.concatenate((sig,zeros))
indices = numpy.tile(numpy.arange(0,frame_len),(numframes,1)) + numpy.tile(numpy.arange(0,numframes*frame_step,frame_step),(frame_len,1)).T
indices = numpy.array(indices,dtype=numpy.int32)
frames = padsignal[indices]
win = numpy.tile(winfunc(frame_len),(numframes,1))
return frames*win
np.fft.fft
<function numpy.fft.fft(a, n=None, axis=-1, norm=None)>
