数字信号处理综合实验
班级： 姓名： 学号：
一、实验目的
1．掌握MATLAB的程序设计方法；
2．掌握数字信号处理的基本理论和基本方法； 3．掌握语音信号的采集与处理方法；
4．掌握用MATLAB设计FIR和IIR数字滤波器的方法； 5．掌握用MATLAB对信号进行分析和处理的方法；
二、实验原理
y=exp(x): 以e为底的指数。
conj(x): 取x的共轭，即改变x的虚部符号。 real(x): 取复数x的实部。
rand(1,N): 生成在0和1之间均匀分布的随机序列，长度为N。 randn(1,N): 生成正态分布(高斯分布)的随机序列，长度为N。 sound(f,fs): 输入参量是音频数据向量、采样频率和转换位数。
X=fft(x,N): 计算序列x的N点FFT。如果x的长度小于N，则在x后面补零；如果x的长度大于N，则对x进行截取；如果不指定参数N，则以x的实际长度作为FFT的点数。
x=ifft(X,N): 计算序列X的N点IFFT。
Y=fftshift(X): 将序列X分成左右两部分并交换位置。
[N, Wn] = buttord(Wp, Ws, Rp, Rs, 's'): 求巴特沃思滤波器的阶数N和3dB边界角频率Wn。Rp和Rs分别为通带波动δ和阻带衰减At，单位均为dB；'s'表示模拟滤波器设计，如无此参数，则表示数字滤波器设计；Wp和Ws分别表示通带边界角频率Ωc和阻带边界角频率Ωr，其值为标量(低通和高通)或双元素向量(带通和带阻)。
[b, a] = butter(N, Wn, 'ftype', 's'): 根据N和Wn求巴特沃思滤波器的系统函数H(s)，b和a分别为H(s)分子和分母多项式的系数。‘ftype’指定滤波器类型，值为‘low’、‘high’或‘stop’， ‘low’为缺省值，表示设计低通或带通滤波器(取决于Wn是标量还是双元素向量)，‘high’或‘stop’分别表示设计高通和带阻滤波器。低通和高通滤波器的阶数为N，带通和带阻滤波器的阶数为2N。
[z, p, k] = butter(N, Wn, 'ftype', 's'): 根据N和Wn求巴特沃思滤波器的零、极点和增益。z、p和k分别为H(s)的零点向量、极点向量和增益。系统函数H(s)与b和a或z、p和k的关系为：
b(1)sM?b(2)sM?1?b(M)s?????b(M?1)H(s)?sN?a(2)sN?1?a(N)s?????a(N?1)k(s?z(1))(s?z(2))?(s?z(M))?(s?p(1))(s?p(2))?(s?p(N))
[N, Wp] = cheb1ord(Wp, Ws, Rp, Rs, 's'): 求切比雪夫I型滤波器的阶数N。输入
参数的含义与巴特沃思滤波器相同，输出参数Wp等于输入参数Wp。
[b, a] = cheby1(N, Rp, Wp, 'ftype', 's'): 求切比雪夫I型滤波器的系统函数H(s)，b和a分别为H(s)分子和分母多项式的系数。
[z, p, k] = cheby1(N, Rp, Wp, 'ftype', 's'): 求切比雪夫I型滤波器的零、极点和增益。
[N, Wp] = ellipord(Wp, Ws, Rp, Rs,'s'): 求椭圆滤波器的阶数N。
[b, a] = ellip(N, Rp, Rs, Wp, 'ftype', 's'): 求椭圆滤波器的系统函数H(s)，b和a分别为H(s)分子和分母多项式的系数。
[z, p, k] = ellip(N, Rp, Rs, Wp, 'ftype', 's'): 求椭圆滤波器的零、极点和增益。 [bz, az] = impinvar(b, a, fs): 通过脉冲响应不变法求数字滤波器的系统函数H(z)。b和a分别为模拟滤波器系统函数H(s)的分子和分母多项式系数；fs为采样频率；bz和az分别为H(z)的分子和分母多项式系数。
[bz, az] = bilinear(b, a, fs): 通过双线性变换法求数字滤波器的系统函数H(z)。输入、输出参数的含义与impinvar函数相同。系统函数H(z)与bz和az的关系为：
bz(1)?bz(2)z?1?bz(3)z?2?????bz(N?1)z?NH(z)?az(1)?az(2)z?1?az(3)z?2?????az(N?1)z?N window=ones(1, N): 产生N点矩形窗，行向量。 window=hann(N): 产生N点汉宁窗，列向量。
window=hanning(N): 产生N点非零汉宁窗，列向量。等价于去除hann(N+2)的第一个零元素和最后一个零元素，得到的N点非零窗函数。
window=hamming(N): 产生N点海明窗，列向量。
window=blackman(N): 产生N点布莱克曼窗，列向量。
window=kaiser(N, beta): 产生参数为beta的N点凯塞窗，列向量。
[M, Wd, beta, ftype]=kaiserord(f, a, dev, fs): 凯塞窗参数估计。f为一组边界频率，最高频率为fs/2。a为f中各个频带的幅度值，通带取1，阻带取0。如果f中有2个元素，则形成3个频带，其中第1个和第3个是通带或阻带，第2个是过渡带，a中也有2个元素，指明第1个和第3个频带是通带还是阻带；如果f中有4个元素，则形成5个频带，其中1，3和5是通带或阻带，2和4是过渡带，a中有3个元素，指明1，3和5是通带还是阻带。dev的维数与a相同，指明每个频带上的波动值。fs为采样频率。M为FIR滤波器的阶数，M=N-1。Wd为归一化边界频率，等于数字边界角频率除以π，或者边界频率除以fs/2。beta就是凯塞窗的参数β。ftype为滤波器的类型。
b = fir1(M, Wd, 'ftype', window): 用窗函数法求FIR滤波器的系数b(单位脉冲响应)。M为滤波器的阶数，M=N-1。Wd为一组归一化边界频率，通带和阻带间隔分布，无过渡带；只有一个元素，表示低通或高通滤波器；有两个元素表示带通和带阻滤波器；有三个及以上元素，表示多带滤波器。'ftype'表示滤波器类型，'high'表示高通滤波器，'stop'表示带阻滤波器，'DC-0'表示多带滤波器的第一个频带为阻带，'DC-1'表示多带滤波器的第一个频带为通带。window为窗口类型，缺省为海明窗。
b = fir2(M, f, m, window): 用频率采样法求FIR滤波器的系数b。M为滤波器的阶数，M=N-1。f为一组归一化频率，第一个元素必须为0，最后一个元素必须为1(对应奈奎斯特频率，即采样频率的一半)，中间的元素按升序排列。m的维数与f相同，指明f中每个频率上的理想幅度。window为窗口类型，缺省为海明窗。Fir2可以实现任意幅度特性的滤波器。
三、实验内容
1、连续信号的频率分析
一个连续信号x(t)含两个频率分量：x(t)?sin(2000πt)?sin(2250πt)，用fs?8000Hz的采样器对该信号进行采样，对采样后的序列分别用长度为16和128的FFT观察其频谱，你观察到什么现象？请解释原因。 2、数字滤波器的设计
分别用双线性变换法和窗口设计法设计IIR和FIR数字滤波器，包括低通、高通和带通三种类型。
三种滤波器的性能指标如下：
（1）低通滤波器：通带和阻带边界频率分别为1000Hz和1200Hz，通带波动为1dB，阻带最小衰减为50dB，采样频率为8000Hz；
（2）高通滤波器：通带和阻带边界频率分别为3000Hz和2800Hz，通带波动为1dB，阻带最小衰减为50dB，采样频率为8000Hz；
（3）带通滤波器：通带边界频率分别为1200Hz和3000Hz，阻带边界频率分别为1000Hz和3200Hz，通带波动为1dB，阻带最小衰减为50dB，采样频率为8000Hz。
请给出每种滤波器的阶数、系统函数和对数幅频响应。 3、语音信号的采集与处理
（1）用Windows自带的录音机或其他录音软件录制一段语音。录制格式如下：采样频率8000Hz，8位，单声道。录音内容为你的姓名和学号，时长控制在5秒左右，说话不要有过多的停顿。
（2）对语音信号进行频谱分析，请画出时域波形和对数幅频特性。
（3）分别用“2、数字滤波器的设计”中的三种FIR滤波器对语音信号进行滤波，比较滤波前后语音信号的时域波形与幅频特性。并用MATLAB的sound函数回放语音，试听滤波前后语音的变化。请在实验结果与分析中解释实验现象。
（4）分别用“2、数字滤波器的设计”中的IIR低通滤波器和线性相位FIR低通滤波器对语音信号进行滤波，比较两种滤波器滤波后的时域波形与幅频特性。并用MATLAB的sound函数回放语音，试听两种滤波器滤波后语音的差异。请在实验结果与分析中解释实验现象。
四、实验结果与分析
1.
程序：
n1=0:1:15;
x=sin(0.25*pi*n1)+sin(0.28125*pi*n1); subplot(2,2,1); plot(n1,x,'-*'); xlabel('n'); ylabel('x1(n)');
title( 'N=16,时域特性'); subplot(2,2,2);
xk1=abs(fft(x)); stem(n1,xk1) xlabel('k'); ylabel('X1(k)');
title( 'N=16,幅频特性'); n2=0:1:127;
x=sin(0.25*pi*n2)+sin(0.28125*pi*n2); subplot(2,2,3); plot(n2,x,'-*'); xlabel('n'); ylabel('x2(n)');
title( 'N=128,时域特性'); subplot(2,2,4); xk2=abs(fft(x)); stem(n2,xk2) xlabel('k'); ylabel('X2(k)');
title( 'N=128,幅频特性'); 结果：
N=16,时域特性210-1-20510151510N=16,幅频特性X1(k)500x1(n)51015nN=128,时域特性210-1-2050n1001508060kN=128,幅频特性X2(k)x2(n)40200050k100150
分析：当采样点N2=128时，幅频特性的拟合性更好。因为所取的N值越大，即采样的点越多，就会越贴近实际，误差也就越小。 2.
(1)低通滤波器：
IIR数字滤波器：
程序： clear;
fs=8000;fc=1000;fr=1200;rp=1;rs=50; wp=2*fs*tan(2*pi*fc/fs/2); ws=2*fs*tan(2*pi*fr/fs/2);
[N, wn]=ellipord(wp, ws, rp, rs, 's'); [b a]=ellip(N,rp,rs,wn,'s'); [bz,az]=bilinear(b,a,fs); [h,w]=freqz(bz,az); f=w/(2*pi)*fs;
plot(f,20*log10(abs(h))) axis([0,3000,-120,10]);
grid; xlabel('频率/Hz'); ylabel('幅度'); legend('双线性变换法'); title('椭圆低通滤波器');
结果：
椭圆低通滤波器 0双线性变换法-20-40幅度-60-80-100-120 050010001500频率/Hz200025003000
系统函数为：
错误！未找到引用源。