用 MATLAB 处理音频 一、 问题的提出： 数字语音是信号的一种，我们处理数字语音信号，也就是对一种信号的处理， 那信号是什么呢？信号是传递信息的函数。离散时间信号——序列——可以用图 形来表示。按信号特点的不同，信号可表示成一个或几个独立变量的函数。例如， 图像信号就是空间位置（二元变量）的亮度函数。一维变量可以是时间，也可以 是其他参量，习惯上将其看成时间。信号有以下几种： （1）连续时间信号：在连续时间范围内定义的信号，但信号的幅值可以是 连续数值，也可以是离散数值。当幅值为连续这一特点情况下又常称为模拟信号。 实际上连续时间信号与模拟信号常常通用，用以说明同一信号。 （2）离时间信号：时间为离散变量的信号，即独立变量时间被量化了。而 幅度仍是连续变化的。 （3）数字信号：时间离散而幅度量化的信号。 语音信号是基于时间轴上的一维数字信号，在这里主要是对语音信号进行 频域上的分析。在信号分析中，频域往往包含了更多的信息。对于频域来说，大 概有 8 种波形可以让我们分析：矩形方波，锯齿波，梯形波，临界阻尼指数脉冲 波形，三角波，余旋波，余旋平方波，高斯波。对于各种波形，我们都可以用一 种方法来分析，就是傅立叶变换：将时域的波形转化到频域来分析。 于是，本课题就从频域的角度对信号进行分析，并通过分析频谱来设计出合 适的滤波器。当然，这些过程的实现都是在 MATLAB 软件上进行的，MATLAB 软件在数字信号处理上发挥了相当大的优势。 二、 设计方案： 利用 MATLAB 中的 wavread 命令来读入（采集）语音信号，将它赋值给某 一向量。再将该向量看作一个普通的信号，对其进行 FFT 变换实现频谱分析， 再依据实际情况对它进行滤波。对于波形图与频谱图（包括滤波前后的对比图） 都可以用 MATLAB 画出。我们还可以通过 sound 命令来对语音信号进行回放， 以便在听觉上来感受声音的变化。 选择设计此方案，是对数字信号处理的一次实践。在数字信号处理的课程学 习过程中，我们过多的是理论学习，几乎没有进行实践方面的运用。这个课题正 好是对数字语音处理的一次有利实践，而且语音处理也可以说是信号处理在实际 应用中很大众化的一方面。 这个方案用到的软件也是在数字信号处理中非常通用的 MATLAB 软件。所 以这个课题的设计过程也是一次数字信号处理在 MATLAB 中应用的学习过程。 课题用到了较多的 MATLAB 语句，而由于课题研究范围所限，真正与数字信号 有关的命令函数却并不多。 三、 主体部分： （一）、语音的录入与打开： [y,fs,bits]=wavread('Blip',[N1 N2]); 用于读取语音，采样值放在向量 y 中，fs 表示采样频率(Hz)，样 bits 表示采 样位数。[N1 N2]表示读取从 N1 点到 N2 点的值（若只有一个 N 的点则表示读 取前 N 点的采值）。 sound(x,fs,bits); 

用于对声音的回放。向量 y 则就代表了一个信号（也即一个复杂的“函数表 达式”）也就是说可以像处理一个信号表达式一样处理这个声音信号。 FFT 的 MATLAB 实现 在 MATLAB 的信号处理工具箱中函数 FFT 和 IFFT 用于快速傅立叶变换和 逆变换。下面介绍这些函数。 函数 FFT 用于序列快速傅立叶变换。 函数的一种调用格式为： y=fft(x) 其中，x 是序列，y 是序列的 FFT，x 可以为一向量或矩阵，若 x 为一向量， y 是 x 的 FFT。且和 x 相同长度。若 x 为一矩阵，则 y 是对矩阵的每一列向量进 行 FFT。 如果 x 长度是 2 的幂次方，函数 fft 执行高速基－2FFT 算法；否则 fft 执行 一种混合基的离散傅立叶变换算法，计算速度较慢。 函数 FFT 的另一种调用格式为 ： y=fft(x,N) 式中，x，y 意义同前，N 为正整数。 函数执行 N 点的 FFT。若 x 为向量且长度小于 N，则函数将 x 补零至长 度 N。若向量 x 的长度大于 N，则函数截短 x 使之长度为 N。若 x 为矩阵，按 相同方法对 x 进行处理。 经函数 fft 求得的序列 y 一般是复序列，通常要求其幅值和相位。MATLAB 提供求复数的幅值和相位函数：abs，angle，这些函数一般和 FFT 同时使用。 函数 abs(x)用于计算复向量 x 的幅值，函数 angle(x)用于计算复向量的相角，介 于 和 之间，以弧度表示。 函数 unwrap(p)用于展开弧度相位角 p ,当相位角绝对变化超过 时，函数 把它扩展至 。 用 MATLAB 工具箱函数 fft 进行频谱分析时需注意： （1） 函数 fft 返回值 y 的数据结构对称性 若已知序列 x=[4,3,2,6,7,8,9,0],求 X(k)=DFT[x(n)]。 利用函数 fft 计算，用 MATLAB 编程如下： N=8; n=0:N-1; xn=[4 3 2 6 7 8 9 0]'; XK=fft(xn) 结果为： XK = 39.0000 -10.7782 + 6.2929i 0 - 5.0000i 4.7782 - 7.7071i 5.0000 4.7782 + 7.7071i 0 + 5.0000i -10.7782 - 6.2929i 

由程序运行所得结果可见，X(k)和 x(n)的维数相同，共有 8 个元素。X(k) 的第一行元素对应频率值为 0，第五行元素对应频率值为 Nyquist 频率，即标准 频率为 1.因此第一行至第五行对应的标准频率为 0～1。而第五行至第八行对应 的是负频率，其 X(k)值是以 Nyquist 频率为轴对称。（注：通常表示为 Nyquist 频率外扩展，标以正值。） 一般而言，对于 N 点的 x(n)序列的 FFT 是 N 点的复数序列，其点 n=N/2+1 对应 Nyquist 频率，作频谱分析时仅取序列 X(k)的前一半，即前 N/2 点即可。X(k)的 后一半序列和前一半序列时对称的。 （2） 频率计算 若 N 点序列 x(n)(n=0,1,…,N-1)是在采样频率 下获得的。它的 FFT 也是 N 点序列，即 X(k)(k=0,1,2,…,N-1)，则第 k 点所对应实际频率值为 f=k*f /N. （3） 作 FFT 分析时，幅值大小与 FFT 选择点数有关，但不影响分析结 果。 2、设计内容： （1）下面的一段程序是语音信号在 MATLAB 中的最简单表现，它实现了语音的 读入打开，以及绘出了语音信号的波形频谱图。 [x,fs,bits]=wavread('ding.wav',[1024 5120]); sound(x,fs,bits); X=fft(x,4096); magX=abs(X); angX=angle(X); subplot(221);plot(x);title('原始信号波形'); subplot(222);plot(X); title('原始信号频谱'); subplot(223);plot(magX);title('原始信号幅值'); subplot(224);plot(angX);title('原始信号相位'); 程序运行可以听到声音，得到的图形为： （2）定点分析：已知一个语音信号，数据采样频率为 100Hz，试分别绘制 N＝ 128 点 DFT 的幅频图和 N＝1024 点 DFT 幅频图。 编程如下： x=wavread('ding.wav'); sound(x); fs=100;N=128; y=fft(x,N); magy=abs(y); f=(0:length(y)-1)'*fs/length(y); subplot(221);plot(f,magy); xlabel('频率（Hz）');ylabel('幅值'); title('N=128(a)');grid subplot(222);plot(f(1:N/2),magy(1:N/2)); xlabel('频率（Hz）');ylabel('幅值'); title('N=128(b)');grid fs=100;N=1024; y=fft(x,N); magy=abs(y); 

f=(0:length(y)-1)'*fs/length(y); subplot(223);plot(f,magy); xlabel('频率（Hz）');ylabel('幅值'); title('N=1024(c)');grid subplot(224);plot(f(1:N/2),magy(1:N/2)); xlabel('频率（Hz）');ylabel('幅值'); title('N=1024(d)');grid （3）若信号长度 T=25.6s，即抽样后 x(n)点数为 T/Ts=256，所得频率分辨率为 Hz，以此观察数据长度 N 的变化对 DTFT 分辨率的影响： 编程如下： [x,fs,bits]=wavread('ding.wav'); N=256; f=0:fs/N:fs/2-1/N; X=fft(x); X=abs(X); subplot(211) plot(f(45:60),X(45:60));grid xlabel('Hz'),ylabel('|H(ejw)|') %数据长度 N 扩大 4 倍后观察信号频谱 N=N*4; f=0:fs/N:fs/2-1/N; X=fft(x); X=abs(X); subplot(212) plot(f(45*4:4*60),X(4*45:4*60));grid xlabel('Hz'),ylabel('|H(ejw)|')