这篇文章上次修改于 371 天前,可能其部分内容已经发生变化,如有疑问可询问作者。

引言

傅里叶变换是熟悉的频域分析手段,在信号分析处理中起到重要的作用。
标准的傅立变换适用于具有稳定周期性变换和平稳随机信号中,但不适合直接分析语音信号。语音信号可以近似地认为在短时间内稳定不变,因此可以分析短时语音研究频域特性。

短时谱分析

取出一段语音信号,进行加窗并取出短时数据帧

close all
clear all
clc
[s,fs]=audioread('D:\englishapp\MATLAB\single\myvoice.wav') ;
s=s(:,1)';
s=endpoint_v(s);
N=512;
h=hamming(N);
ss=s(5*N+1:6*N).*h'; 
sf=20*log(abs(fft(ss)));
subplot(2,1,1); 
plot([1:N]/fs,ss);
title('original signal');
grid on
subplot(2,1,2)
plot(linspace(1,fs/2,N/2),sf(1:N/2));
title('短时谱');
grid on

运行结果:
[RiLwr9.png]

分析:如图所示,在窗长为N=256的情况下,从语音的短时谱上可以看出,短时频谱图中存在一定程度的慢变化,这些是声道滤波器的共振峰,从一个基音周期到下一个基音周期,共振峰会发生变化,如2500hz到3500hz的共振峰与1500hz到2500hz有显著的不同。
在短时谱的截止频率处比较尖锐,这是由于汉明窗具有低通的性质。

语谱图分析

close all
clear all
clc
[s,fs]=audioread('D:\englishapp\MATLAB\single\myvoice.wav') ;
s=s(:,1)';
s=endpoint_v(s);
specgram(s,512,fs,100);  %specgram(y,nfft,Fs,window,overlap)
xlabel('时间(s)'); 
ylabel('频率(Hz)'); 
title('语谱图');

运行结果:
[RiL0bR.png]
分析:从语谱图上看,在0.1s到0.45s,2000hz到4000hz之间有较多的能量集中,而0-500hz的能量几乎占据持续了全程。这些能量集中的地方在整个语音中所占的比重较大,集中的能量在时间轴中形成横向的纹路,这些纹路就是声纹。声纹具有个体差异性,在语音信号的识别、分析处理、司法鉴别、安全等等方面有广泛的作用。

倒谱和复倒谱分析

% claculation complex cepstrum with the Differentiation
 
close all;  clear all;  clc
 
 
N = 256 * 1;    INC = 128;  nfft = N * 1;
[s fs] = audioread('D:\englishapp\MATLAB\single\myvoice.wav');
s=s(:,1)';  ss = enframe(s, N, INC);   
ns = ss(53,:);  nns = [1: length(ns)] .* ns;
fns = fft(ns, nfft);    fnns = fft(nns, nfft);
 
fth = fnns ./ fns;  ifth = ifft(fth);
hat_ns = real(ifth) ./ [1: length(ns)];
 
h=hamming(N)'
%figure; plot(hat_ns)
figure; 
subplot(2,2,1),plot(cceps(ns)),title('复倒谱')
subplot(2,2,2),plot(rceps(ns)),title('倒谱')
subplot(2,2,3),plot(cceps(ns).*h),title('N=256汉明窗复倒谱')
subplot(2,2,4),plot(rceps(ns).*h),title('N=256汉明窗倒谱')

[RFiePJ.png]

分析:如图为使用微分法计算的倒谱和复倒谱,在加入N=256的汉明窗之后,倒谱和复倒谱的共振峰更加的明显,虽然幅度降低了,但是各个峰可以清晰地看出来。以此来进行基音检测和共振峰检测更加容易。
在基音检测上,倒谱的高频部分包含了基音周期的峰值信息,如图所示,倒谱的高频部分为x=89,因此,基因周期应该为89hz。
在共振峰检测上,平滑化过的谱的峰值对应共振峰的频率。以此可以估计共振峰。

附录

endpoint_v.m文件:

function [ output_args ] = endpoint_v( input_args )
%通过端点检测,提取有效的数据段
x=input_args;
%幅度归一化到[-1,1]
% x = double(x);
x = x / max(abs(x));

%常数设置
FrameLen = 256;
FrameInc = 64;

amp1 = 10;
amp2 = 2;
% amp2=0.5;
zcr1 = 10;
zcr2 = 5;

maxsilence = 8;  
minlen  = 15;    
status  = 0;
count   = 0;
silence = 0;

%计算过零率
tmp1  = enframe(x(1:end-1), FrameLen, FrameInc);
tmp2  = enframe(x(2:end)  , FrameLen, FrameInc);
signs = (tmp1.*tmp2)<0;
diffs = (tmp1 -tmp2)>0.02;
zcr   = sum(signs.*diffs, 2);

%计算短时能量
amp = sum(abs(enframe(filter([1 -0.9375], 1, x), FrameLen, FrameInc)), 2);

%调整能量门限
amp1 = min(amp1, max(amp)/4);
amp2 = min(amp2, max(amp)/8);

%开始端点检测
x1 = 0; 
x2 = 0;
for n=1:length(zcr)
   switch status
   case {0,1}                        % 0 = 静音, 1 = 可能开始
      if amp(n) > amp1          % 确信进入语音段
         x1 = max(n-count-1,1);
         status  = 2;
         silence = 0;
         count   = count + 1;
      elseif amp(n) > amp2   % 可能处于语音段
          status = 1;
          count  = count + 1;
      else                                 % 静音状态
         status  = 0;
         count   = 0;
      end
   case 2,                               % 2 = 语音段
      if amp(n) > amp2           % 保持在语音段
         count = count + 1;
      else                                % 语音将结束
         silence = silence+1;
         if silence < maxsilence % 静音还不够长,尚未结束
            count  = count + 1;
         elseif count < minlen   % 语音长度太短,认为是噪声
            status  = 0;
            silence = 0;
            count   = 0;
         else                               % 语音结束
            status  = 3;
         end
      end
   case 3,
      break;
   end
end

count = count-silence/2;
x2 = x1 + count -1;
output_args=input_args(x1*FrameInc:x2*FrameInc);