河海大学数字信号处理实验 综合实验

数字信号处理综合实验

班级： 姓名： 学号：

一、实验目的

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

系统函数为：

错误！未找到引用源。

FIR数字滤波器：

通带和阻带的数字边界频率:

?c?0.25??r?0.3?

求理想低通滤波器的边界频率用通带和阻带边界频率的中心近似:

?n??c??r2?0.275?

滤波器的过渡带宽为 ?r??c?0.05?， 阻带衰减不小于50dB，因此选择海明窗，因此窗口长度为:

6.6??0.05??N?132 N线性相位延迟常数为:

??N?1?65.5 2根据理想边界频率?n和线性相位延迟常数?，求理想单位脉冲响应hd(n):

1hd(n)?2?????nej?(n??)d??nsin[(n??)?n]

?(n??)窗函数与理想单位脉冲响应相乘，得到线性相位FIR低通滤波器的单位脉冲响应:

????2?n??sin??n?65.5?0.275hd(n)??0.54?0.46cos??R132?n? ????131?n?65.5????程序：

N=132; n=0:N-1; a=(N-1)/2; wn=0.275*pi;

hd=sin(0.275*pi*(n-a))./(pi*(n-a));

win=0.54-0.46*cos(n*2*pi/(N-1)); %win=hamming(N)'; h=win.*hd;

figure; stem(n,h);

xlabel('n'); ylabel('h(n)'); grid;

title('线性相位FIR低通滤波器的单位脉冲响应h(n)'); [H,w]=freqz(h,1); H=20*log10(abs(H)); figure; plot(w/pi, H);

axis([0 1 -100 10]); xlabel('频率/Hz'); ylabel('幅度/dB'); grid;

title('线性相位FIR低通滤波器，海明窗，N=132');

结果：

线性相位FIR低通滤波器的单位脉冲响应h(n)0.30.250.20.15h(n)0.10.050-0.05-0.10204060n80100120140

线性相位FIR低通滤波器，海明窗，N=132100-10-20-30幅度/dB-40-50-60-70-80-90-10000.10.20.30.40.50.6频率/Hz0.70.80.91

(2)高通滤波器：

IIR数字滤波器：

程序： clear

fs=8000;fr=2800;fc=3000;rs=50;rp=1; wp=2*fs*tan(2*pi*fc/fs/2); ws=2*fs*tan(2*pi*fr/fs/2);

[N, wp] = ellipord(wp, ws, rp, rs, 's'); [b,a]=ellip(N,rp,rs,wp,'high','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,4000,-100,10]);grid; xlabel('频率/Hz'); ylabel('幅度'); legend('双线性变换法'); title('椭圆高通滤波器');

结果：

椭圆高通滤波器100-10-20-30双线性变换法 幅度-40-50-60-70-80-90-100 0500100015002000频率/Hz2500300035004000

系统函数：

错误！未找到引用源。

N=6

FIR滤波器：

通带和阻带的数字边界频率:

?c?0.75??r?0.7?

理想高通滤波器的边界频率用通带和阻带边界频率的中心近似:

?n??c??r2?0.725?

滤波器的过渡带宽为 ?c??r?0.05?， 阻带衰减不小于50dB，因此选择海明窗因此窗口长度为:

6.6??0.05??N?132 N线性相位延迟常数为:

??N?1?65.5 2根据理想边界频率?n和线性相位延迟常数?，求理想单位脉冲响应hd(n):

hd(n)?12??????nej?(n??)d??12????nej?(n??)d??sin[(n??)?]?sin[(n??)?n]

?(n??)窗函数与理想单位脉冲响应相乘，得到线性相位FIR低通滤波器的单位脉冲响应:

????2?n??sin??n?65.5????sin??n?65.5?0.725hd(n)??0.54?0.46cos??R132?n? ????n?65.5??131???程序：

N=132; n=0:N-1; a=(N-1)/2; wn=0.725*pi;

hd=(sin(pi*(n-a))-sin(0.725*pi*(n-a)))./(pi*(n-a)); win=0.54-0.46*cos(n*2*pi/(N-1)); %win=hamming(N)'; h=win.*hd;

figure; stem(n,h);

xlabel('n'); ylabel('h(n)'); grid;

title('线性相位FIR高通滤波器的单位脉冲响应h(n)'); [H,w]=freqz(h,1); H=20*log10(abs(H)); figure; plot(w/pi, H);

axis([0 1 -100 10]); xlabel('频率/Hz'); ylabel('幅度/dB'); grid; title('线性相位FIR高通滤波器，海明窗，N=132');

结果：

线性相位FIR高通滤波器的单位脉冲响应h(n)0.20.150.10.05h(n)0-0.05-0.1-0.15-0.20204060n80100120140线性相位FIR高通滤波器，海明窗，N=132100-10-20-30幅度/dB-40-50-60-70-80-90-10000.10.20.30.40.50.6频率/Hz0.70.80.91(3)带通滤波器 IIR：

程序： clear;

fs=8000;fc=[1200 3000];fr=[1000 3200];rp=1;rs=50; wp=2*fs*tan(2*pi*fc/fs/2); ws=2*fs*tan(2*pi*fr/fs/2);

[N, wn] = cheb1ord(wp, ws, rp, rs, 's'); [b2 a2]=cheby1(N,rp,wp,'s'); [bz2,az2]=bilinear(b2,a2,fs); [h2,w]=freqz(bz2,az2); f=w/(2*pi)*fs;

figure; plot(f,20*log10(abs(h2)),'-b'); grid; xlabel('频率/Hz'); ylabel('幅度'); legend('双线性变换法');

title('切比雪夫带通滤波器');

切比雪夫带通滤波器0双线性变换法-50 -100幅度-150-200-250-300 0500100015002000频率/Hz2500300035004000

系统函数：

H(z)?1.0000-1.3195z?1?4.5855z?2?5.1873z?3?11.4866z?4-11.3175z?5?19.0861z?6?0.0092z?7?0.1368z?8?0.0104z?91?0.8124z?1?2.1656z?2?1.6670z?3?3.6729z?4?2.2977z?5?3.5841z?6?1.8949z?7?2.5650z?8?0.9933z?9阶数N=9

FIR：

求理想通带和阻带的数字边界频率和过渡带宽:

?c1?2??r1?2??1?fc1?0.3?fsfr1?0.25?fs2?0.275??c2?2??r2?2??2?fc2?0.75?fsfr2?0.8?fs2?0.775?

?c1??r1?c2??r2滤波器的过渡带宽为 ?r1??c1?0.05?

阻带衰减不小于50dB，因此选择汉明窗，因此窗口长度为:

6.6??0.05??N?132 N线性相位延迟常数为:

??N-1?65.5 2根据理想边界频率?n和线性相位延迟常数?，求理想单位脉冲响应hd(n):

hd(n)?12??????12ej?(n??)d??12????21ej?(n??)d??sin[(n??)?2]?sin[(n??)?1]

?(n??)窗函数与理想单位脉冲响应相乘，得到线性相位FIR低通滤波器的单位脉冲响应:

????sin??n?65.5?0.275???2?n??sin??n?65.5?0.775hd(n)??0.54?0.46cos??R132?n? ????n?65.5??131???程序：

N=132; n=0:N-1;

a=(N-1)/2; fs=20000; w1=0.275*pi; w2=0.775*pi;

hd=(sin(w2*(n-a))-sin(w1*(n-a)))./(pi*(n-a));

win=0.54-0.46*cos(n*2*pi/(N-1)); %window=hamming(N); h=win.*hd;

figure; stem(n,h);

xlabel('n'); ylabel('h(n)'); grid;

title('线性相位FIR带通滤波器的单位脉冲响应h(n)'); [H,w]=freqz(h,1); H=20*log10(abs(H)); figure; plot(w/(2*pi)*fs, H);

axis([0 0.5*fs -100 10]); xlabel('频率/Hz'); ylabel('幅度/dB'); grid;

title('线性相位FIR带通滤波器，海明窗，N=132'); 结果：

线性相位FIR带通滤波器的单位脉冲响应h(n)0.40.30.20.1h(n)0-0.1-0.2-0.3-0.40204060n80100120140线性相位FIR带通滤波器，海明窗，N=132100-10-20-30幅度/dB-40-50-60-70-80-90-1000100020003000400050006000频率/Hz700080009000100003.

（2）

程序：

[x fs bits]=wavread('F:/姓名学号.wav') figure(1) plot(x)

xlabel('t');ylabel('x(t)') title('语音信号时域图像') X=fft(x) figure(2)

stem(20*log10(abs(X)))

xlabel('频率/Hz');ylabel('幅度/dB'); title('语音信号对数频谱')

结果：

语音信号时域图像10.80.60.40.2x(t)0-0.2-0.4-0.6-0.800.511.5t22.533.5x 104

语音信号对数频谱604020幅度/dB0-20-40-6000.511.52频率/Hz2.533.5x 104

（3）

FIR低通滤波： 程序：

[x fs bits]=wavread('F: /姓名学号/.wav') X=fft(x)

N=132;n=0:N-1; a=(N-1)/2; wn=0.275*pi;

hd=sin(wn*(n-a))./(pi*(n-a)); win=0.54-0.46*cos(n*2*pi/(N-1)); h=win.*hd; y=fftfilt(h,x) Y=fft(y) figure(1) subplot(211)

plot(x);title('原始信号时域') subplot(212)

plot(y);title('低通滤波后时域') figure(2) subplot(211)

stem(20*log10(abs(X))) title('原始信号对数幅频') subplot(212)

stem(20*log10(abs(Y)))

本文来源：https://www.bwwdw.com/article/9ph5.html