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

1、数字信号处理综合实验 班级： 姓名： 学号：一、实验目的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): 输入参量是音频数据向量、采样频率和转换

2、位数。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，其值为标量(低通和高通

3、)或双元素向量(带通和带阻)。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)的零点向量、极点向量和增益。系统

4、函数H(s)与b和a或z、p和k的关系为：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 = elli

5、p(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)。输入、输出参数

6、的含义与impinvar函数相同。系统函数H(z)与bz和az的关系为：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=kai

7、serord(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就是凯塞窗

8、的参数。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): 用频率采样法求

9、FIR滤波器的系数b。M为滤波器的阶数，M=N-1。f为一组归一化频率，第一个元素必须为0，最后一个元素必须为1(对应奈奎斯特频率，即采样频率的一半)，中间的元素按升序排列。m的维数与f相同，指明f中每个频率上的理想幅度。window为窗口类型，缺省为海明窗。Fir2可以实现任意幅度特性的滤波器。三、实验内容1、连续信号的频率分析一个连续信号含两个频率分量：，用Hz的采样器对该信号进行采样，对采样后的序列分别用长度为16和128的FFT观察其频谱，你观察到什么现象？请解释原因。2、数字滤波器的设计 分别用双线性变换法和窗口设计法设计IIR和FIR数字滤波器，包括低通、高通和带通三种类型。三种滤

10、波器的性能指标如下：（1）低通滤波器：通带和阻带边界频率分别为1000Hz和1200Hz，通带波动为1dB，阻带最小衰减为50dB，采样频率为8000Hz；（2）高通滤波器：通带和阻带边界频率分别为3000Hz和2800Hz，通带波动为1dB，阻带最小衰减为50dB，采样频率为8000Hz；（3）带通滤波器：通带边界频率分别为1200Hz和3000Hz，阻带边界频率分别为1000Hz和3200Hz，通带波动为1dB，阻带最小衰减为50dB，采样频率为8000Hz。请给出每种滤波器的阶数、系统函数和对数幅频响应。3、语音信号的采集与处理（1）用Windows自带的录音机或其他录音软件录制一段语音

11、。录制格式如下：采样频率8000Hz，8位，单声道。录音内容为你的姓名和学号，时长控制在5秒左右，说话不要有过多的停顿。（2）对语音信号进行频谱分析，请画出时域波形和对数幅频特性。（3）分别用“2、数字滤波器的设计”中的三种FIR滤波器对语音信号进行滤波，比较滤波前后语音信号的时域波形与幅频特性。并用MATLAB的sound函数回放语音，试听滤波前后语音的变化。请在实验结果与分析中解释实验现象。（4）分别用“2、数字滤波器的设计”中的IIR低通滤波器和线性相位FIR低通滤波器对语音信号进行滤波，比较两种滤波器滤波后的时域波形与幅频特性。并用MATLAB的sound函数回放语音，试听两种滤波器滤

12、波后语音的差异。请在实验结果与分析中解释实验现象。4、 实验结果与分析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,

13、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,幅频特性);结果：分析：当采样点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

14、*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(椭圆低通滤波器);结果：系统函数为： FIR数字滤波器： 通带和阻带的数字边界频率:求理想低通滤波器的边界频率用通带和阻带边界频率的中心近似

15、:滤波器的过渡带宽为 ，阻带衰减不小于50dB，因此选择海明窗，因此窗口长度为:线性相位延迟常数为:根据理想边界频率和线性相位延迟常数，求理想单位脉冲响应:窗函数与理想单位脉冲响应相乘，得到线性相位FIR低通滤波器的单位脉冲响应:程序：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;titl

16、e(线性相位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);结果：(2) 高通滤波器：IIR数字滤波器：程序：clearfs=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

17、(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(椭圆高通滤波器);结果：系统函数： N=6FIR滤波器： 通带和阻带的数字边界频率:理想高通滤波器的边界频率用通带和阻带边界频率的中心近似:滤波器的过渡带宽为 ，阻带衰减不小于50dB，因此选择海明窗因

18、此窗口长度为:线性相位延迟常数为:根据理想边界频率和线性相位延迟常数，求理想单位脉冲响应:窗函数与理想单位脉冲响应相乘，得到线性相位FIR低通滤波器的单位脉冲响应:程序：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高通滤波器的单位脉冲

19、响应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);结果：(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

20、, 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(切比雪夫带通滤波器);系统函数：阶数N=9FIR： 求理想通带和阻带的数字边界频率和过渡带宽:滤波器的过渡带宽为 阻带衰减不小于50dB，因此选择汉明窗，因此窗口长度为:线性相位延迟常数为:根据理想边界频率和线性相位延迟常数，求理想单位脉冲响

21、应:窗函数与理想单位脉冲响应相乘，得到线性相位FIR低通滤波器的单位脉冲响应:程序：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=

22、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);结果：3.（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(语

23、音信号对数频谱)结果：（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*lo

24、g10(abs(X)title(原始信号对数幅频)subplot(212)stem(20*log10(abs(Y)title(低通滤波后对数幅频)sound(y)结果：FIR高通滤波：程序：x fs bits=wavread(F:/ /姓名学号.wav)X=fft(x)N=132;n=0:N-1;a=(N-1)/2; wn=0.725*pi; hd=sin(pi*(n-a)./(pi*(n-a)-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)subp

25、lot(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)title(高通滤波后对数幅频)sound(y)结果：FIR带通滤波：程序:x fs bits=wavread(F: /姓名学号/.wav)X=fft(x)N=132;n=0:N-1;hd=(sin(0.775*pi*(n-65.5)-sin(0.275*pi*(n-65.5)./(pi*(n-65.5); window=hamming(N)h=window.*hd;y=fftfilt(h,x)Y=fft(y)figure(1)subplot(211)plot(x);title(原始信号时域)subplot(212)plot(y);title(带通滤波后时域)figure(2)subplot(211