信号处理中简单实用的方法——提取信号中的包络

一、用希尔伯特变换计算信号的包络

在求某一信号包络时用得最多的是希尔伯特变换，但并不是希尔伯特变换适用于所有信号求包络的情况。这是因为对于包络没有一个很严格的定义，在求包络时不同的情况会有不同的要求。下面将介绍用希尔伯特变换求取信号的包络。

对MATLAB中自带的希尔伯特变换的函数介绍如下。
名称：hilbert
功能：把序列x(n)作希尔伯特变换为y(n)，又把x(n)和y(n)构成解析信号的序列z(n)=x(n)+jy(n)

调用格式：
z=hilbert(x)
说明：函数hilbert不是单纯地把x(n)作希尔伯特变换得到y(n)，而是得到y(n)后与x(n)
共同构成解析信号序列z(n)，并可以对z(n)直接求模值和相角。
案例1、设信号x(n)=120+96e^[(n/1500)^2]*cos(2π*n/600)，n=-5000：20：5000，求该信号的包络线。设置信号后直接调用hilbert函数求信号的包络，运行第一部分程序如下：

clear all; clc; close all;
 
n=-5000:20:5000;            % 样点设置
% 程序第一部分:直接做做希尔伯特变换
N=length(n);                % 信号样点数
nt=0:N-1;                   % 设置样点序列号
x=120+96*exp(-(n/1500).^2).*cos(2*pi*n/600); % 设置信号
Hx=hilbert(x);              % 希尔伯特变换
% 作图
plot(nt,x,'k',nt,abs(Hx),'r');
grid; legend('信号','包络');
xlabel('样点'); ylabel('幅值')
title('信号和包络')
set(gcf,'color','w');
pause

 运行结果如下：

 

运行程序第一部分后得上图,从图中可以看出通过hilbert函数求信号所得的包络似乎不是信号的包络，而是原有的信号，这是为什么呢？之所以会得到这样一个不理想的包络线完全是直流分量造成的，故对原始信号消除直流分量后再通过hilbert函数求信号的包络。
第二部分程序如下：

% 程序第二部分:消除直流后做希尔伯特变换
y=x-120;                    % 消除直流分量
Hy=hilbert(y);              % 希尔伯特变换
% 作图
figure(2)
plot(nt,y,'k',nt,abs(Hy),'r');
grid; legend('信号','包络');
xlabel('样点'); ylabel('幅值')
title('信号和包络')
set(gcf,'color','w');
figure(3);
plot(nt,x,'k',nt,abs(Hy)+120,'r');
grid; legend('信号','包络'); hold on;
xlabel('样点'); ylabel('幅值')
title('信号和包络')
set(gcf,'color','w');
pause

 运行结果如下：

同时把包络线加上直流分量后叠加在原始信号上，修改了上图的不正确包络线，得下图：

这说明如果信号有直流分量，可以先消除直流分量求出信号的包络线，再在包络线上叠加直流分量恢复原始信号的包络线。

二、用倒谱法来计算语音信号频谱的包络

语音信号频谱的包络线对语音分析来说是比较重要的，它反映了人类发声器官的共振结构，从频谱的包络可提取共振峰参数（频率和带宽）。在语音分析中常用倒谱方法来提取语音信号频谱的包络线。

倒谱（Cepstrum）分析的实质就是对幅值谱取对数后再做一次频谱分析，所以又称为二次频谱分析。

我们知道，若两个信号在频域有不同的频带，则可以通过滤波把这两个信号分离开。同样在倒谱域中，当两个信号在倒谱域中占有不同的倒频率时，也可以通过倒滤波（Lifter）把这两个信号在倒频域上分离。
声带产生的脉冲序列振动频率较高，而口腔是通过肌肉运动发声的，振动频率较低，所以二者在倒谱中处于不同的倒频带内。通过倒滤波把这两个信号分离，可以分别得到激励脉冲在频域的响应和声道在频域的响应。一般把声道在频域的响应称为语音信号频谱的包络线，通过该包络线可提取语音共振峰的信息。

案例2、 从文件su1.txt读入语音数据，求取该语音信号的频谱包络线。程序清单如下：

clear all; clc; close all;
y=load('su1.txt');                            % 读入数据
fs=16000; nfft=1024;                          % 采样频率和FFT的长度
time=(0:nfft-1)/fs;                           % 时间刻度
nn=1:nfft/2; ff=(nn-1)*fs/nfft;               % 计算频率刻度
Y=log(abs(fft(y)));                           % 取幅值的对数
z=ifft(Y);                                    % 按式(4-3-16)求取倒谱
mcep=29;                                      % 分离声门激励脉冲和声道冲击响应
zy=z(1:mcep+1);
zy=[zy' zeros(1,1024-2*mcep-1) conj(zy(end:-1:2))']; % 构建声道冲击响应的倒谱序列
ZY=fft(zy);                                   % 计算声道冲击响应的频谱
% 作图
plot(ff,Y(nn),'k'); hold on;                  % 画出信号的频谱图
plot(ff,real(ZY(nn)),'k','linewidth',2.5);    % 画出包络线
grid; hold off; ylim([-4 5]);
title('信号频谱和声道冲击响频谱（频谱包络）')
ylabel('幅值'); xlabel('频率/Hz'); 
legend('信号频谱','频谱包络')
set(gcf,'color','w');

运行结果如下：

 

实验数据su1.txt下载链接如下：

https://mp.csdn.net/mp_download/manage/download/UpDetailed

参考文献：MATLAB数字信号处理85个实用案例精讲——入门到进阶；宋知用（编著）