语音信号的短时频谱分析-matlab

一、原理

1、语音信号特征

 (1) 在频域内,语音信号的频谱分量主要集中在300-3400Hz的范围内。
利用这个特点,可以用一个防混叠的带通滤波器将此范围内的语音信号频率分出,然后按8kHz的采样
 (2) 在时域内,语音信号具有“短时性”的特点
即在总体上,语音信号的特征是随着时间而变化的，但在一段较短时间间隔内，语音信号保持平稳。在浊音段表现出周期信号的特征，在轻音端表现出随机噪声的特征。

2、语音信号采集的基本原理

        将模拟信号变为数字信号，必须经过采样和量化，得到时间和幅度上均为离散的数字信号语音。采样时， 采样频率必须以高于受测信号的最高频率两倍以上的速度进行取样,才能正确地重建信号。

3、语音分帧

        语音信号是一个准稳态的信号，若把它分成较短的帧，每帧中可将其看做稳态信号，可用处理稳态信号的方法来处理。为了使一帧与另一帧之间的参数能够平稳过渡，应在相邻两帧之间互相有部分重叠。一般情况下，帧长取10 ~ 30ms，所以每秒的帧数约为33 ~ 100帧。帧移与帧长的比值一般取0~1/2。

 设每帧长度为wlen，后一帧对前一帧的位移量为inc，则重叠部分overlap=wlen-inc
MATLAB中的分帧函数：

其中，x为语音信号，win是帧长，inc是帧移。
enframe函数的返回值是一个[帧数 x 帧长]的矩阵，帧数向上取整。

4、窗函数

窗函数是在时域上对信号进行加权的一种处理方法，常用的三种窗函数分别是矩形窗、海宁窗和汉明窗。这三种窗函数的主瓣宽度和第一旁瓣的衰减关系如下：

窗函数	主瓣宽度与频率分辨率的比值	第一旁瓣的衰减（dB）
矩形窗	1	-13
海宁窗	1.54	-31
汉明窗	1.36	-44

从表中可以看出，矩形窗主瓣的宽度 最窄，但第一旁瓣的衰落最小。也就是说，它的频谱泄露要比另外两种窗函数大。在语音分析中，可根据不同情况选择不同的窗函数。

5、短时傅里叶变换

短时傅里叶变换（Short-Time Fourier Transform，STFT）是傅里叶变换的一种扩展，用于对信号的时变特性进行分析。将信号分割为一段一段的小时间窗口，并对每个窗口做傅里叶变换，以得到每个时间段的频率成分。这样可以在时间和频率两个维度上对信号进行分析。STFT通常用于语音信号分析、音频信号处理、图像处理等领域。

短时傅里叶变换 的定义公式如下：

其中，X(n,ω)是在时间窗口n和频率ω上的STFT值；x(m)是输入信号的样本值；w(m−n)是窗口函数，通常使用汉宁窗或矩形窗；N是窗口的长度。

二、设计方案

1、读取音频文件。

   调用audioread，处理得到一个保存音频数据的数组，和一个采样频率

2、确定相关参数

    如窗函数、窗长、重叠点数，重叠长度，傅里叶点数 等
    窗函数:选择海宁窗       
    窗长：win_sz=128
    重叠长度：nooverlap = win_sz - 1;（随着noverlap参数的引入，增大了时间轴分辨率，即每隔（Nw - noverlap）长度进行一次频率轴的更新，随着noverlap逐渐接近Nw，图像上表现为时间轴更加“细腻”，但随之而来的肯定是计算次数的增加。
傅里叶点数：通常取最接近信号长度的2的整数次幂，

     
3、 调用spectrogram函数

做短时傅里叶变换。（S-将输入信号做STFT处理后得到的二维含时间、频率序列的数组数据；更为详细的处理过程见spectrogram函数定义）

4、根据处理后的时频矩阵，绘制语谱图

三、实现代码

Fs = 8000; %设置采样率
R=audiorecorder(Fs,16,1); %录音
recordblocking(R,2); %录制2秒
stop(R)
disp('结束录音');
r = getaudiodata(R); 
rmax=max(abs(r)); %归一化
r=r/rmax;
audiowrite('xinhai.wav',r,Fs); %存储录音
play(R) %播放音频
[Y,Fs]=audioread('xinhai.wav') %读取音频信号
N = length(Y); %求取抽样点数
t= (0:N-1)/Fs; %显示实际时间
y= fft(Y); %对信号进行傅里叶变换
f = Fs/N*(0:round(N/2) - 1); %显示实际频点的一半，频域映射，转化为Hz
subplot(2,1,1);
plot(t,Y);
title('信号时域波形图'); %语音信号的时域图形
xlabel('时间/s');
ylabel('幅度');
subplot(2,1,2)
plot(f,abs(y(1:round(N/2)))); %绘制频域波形
title("信号幅度谱")
xlabel("频率/HZ");
ylabel('幅度');
figure; % 绘图
win_sz = 128; % 窗函数长度设置为128
han_win = hanning(win_sz); % 选择海宁窗
nfft = 2^nextpow2(length(han_win)); % 重叠点数
nooverlap = win_sz - 1; % 重叠长度
[S, F, T, P] = spectrogram(Y(:,1), window, nooverlap, nfft, Fs); 
imagesc(T, F, log10(abs(S))) % 绘图函数
colorbar;
set(gca, 'YDir', 'normal') % 坐标轴设置
xlabel('Time (secs)')
ylabel('Freq (Hz)')
title('STFT transform spectrum')

四、仿真结果

时域图和幅度谱：

语谱图：

 窗设计器：