数字信号处理实验报告

实验一：用FFT做谱分析 1. 实验目的

（1）进一步加深 DFT 算法原理和基本性质的理解。 （2）熟悉FFT算法原理和FFT子程序的应用。

（3）学习用 FFT对连续信号和时域离散信号进行谱分析的方法，了解可能出现的分析误差及其原因，以便在实际中正确应用FFT。

2. 实验环境

MATLAB R2016a。

3. 实验报告要求

（1）简述实验原理及目的。

（2）结合实验中所得给定典型序列幅频特性曲线，与理论结果比较，并分析说明误差产生的原因以及用 FFT 作谱分析时有关参数的选择方法。 （3）总结实验所得主要结论。

4. 实验原理

用FFT对信号做频谱分析是学习数字信号处理的重要内容。经常需要进行谱分析的信号是模拟信号和时域离散信号。对信号进行谱分析的重要问题是频谱分辨率D和分析误差。频谱分辨率直接和FFT的变换区间N有关，因为FFT能够实

π/ND。可根据此式选择FFT的变换区间N。误差主要来现的频率分辨率是2自于用FFT作频谱分析时，得到的是离散谱，而信号（周期信号除外）是连续谱，只有当N较大时，离散谱的包络才能逼近于连续谱，因此N要适当选择大一些。

周期信号的频谱是离散谱，只有用整数倍周期的长度作FFT，得到的离散谱才能代表周期信号的频谱。如果不知道信号周期，可以尽量选择信号的观察时间长一些。

对模拟信号进行谱分析时，首先要按照采样定理将其变成时域离散信号。如果是模拟周期信号，也应该选取整数倍周期的长度，经过采样后形成周期序列，按照周期序列的谱分析进行。

5. 实验内容

（1）编制信号产生子程序，产生以下典型信号供谱分析用。

第 1 页

x1(n)R4(n)n10n3 x2(n)8n4n7 x5(n)sinn

80otherwisex6(t)cos8tcos16tcos20t4n0n3x3(n)n34n70otherwise应该注意，如果给出的信号是连续信号时（如x6(t)），则首先要根据其最高频率确定采样频率Fs以及由频率分辨率选择采样点数 N，然后对其进行软件采样，产生对应序列x(n)x(nT),0nN1。其中T1/Fs。对于信号x6(t)，频率分辨率的选择要以能分辨开来其中的三个频率对应的谱线为准则。对周期序列，最好截取周期的整数倍进行谱分析，否则有可能产生较大的分析误差。 （2）编写主程序，对（1）中所给出的信号逐个进行谱分析。 （3）按实验内容要求，上机实验，并写出实验报告。 实验代码： clc,clear;

x1n=ones(1,4); %x1=R4(n)n1=0:3; subplot 231;stem(n1,x1n,'.');

xlabel('n');ylabel('x1(n)');title('x1(n)=R4(n)');axis([-1,4,0,2]); %============================================================== n2=0:7;x2n=[1,2,3,4,4,3,2,1]; subplot 232;stem(n2,x2n,'.');

xlabel('n');ylabel('x2(n)');title('x2(n)');axis([-1,8,0,5]); %============================================================== n3=0:7;x3n=[4,3,2,1,1,2,3,4]; subplot 233;stem(n3,x3n,'.');

xlabel('n');ylabel('x3(n)');title('x3(n)');axis([-1,8,0,5]); %============================================================== n4=0:7;x4n=cos(n4*pi/4); subplot 234;stem(n4,x4n,'.');

xlabel('n');ylabel('x4(n)');title('x4(n)');axis([-1,8,-2,2]); %============================================================== n5=0:15;x5n=sin(pi/8*n5); subplot 235;stem(n5,x5n,'.');

xlabel('n');ylabel('x5(n)');title('x5(n)');axis([-1,16,-2,2]); %============================================================= n6=0:0.5/16:0.5;

x6n=cos(8*pi*n6)+cos(16*pi*n6)+cos(20*pi*n6); subplot 236;

stem(n6,x6n,'.');xlabel('n');ylabel('x6(n)');title('x6(n)');

第 2 页

x4(n)cosn4

实验结果：

对于时域连续信号x6(t)cos8tcos16tcos20t,取采样频率

Fs32Hz20Hz，采样点N=16，分辨率FFs/N2Hz，TpN/Fs0.5s为一个周期。

谱分析实验代码：

编写DFT函数myfft如下：

function Xk=myfft(xn,N)

%编写myfft实现DFT进行频谱分析

k=0:N-1;xn=[xn,zeros(1,N-length(xn))]; Xk=zeros(1,N);a=0:N-1; for n=1:N for m=1:N

Xk(n)=Xk(n)+xn(m)*exp(-j*2*pi/N*k(n)*a(m)); end end

分别作出各个信号的8点DFT频谱图线与FFT变换后的理想幅频响应图线 （程序以x1(n)为例）： clc,clear;

k1=0:7;k2=0:15;

%%%%%%%%%%%%%%%%%%%%%%%%x1(n)%%%%%%%%%%%%%%%%%%%%%%%%%% x1n=ones(1,4);n1=0:3;

subplot 311;stem(n1,x1n,'.');axis([-1,4,0,1.5]); xlabel('n');ylabel('x1(n)');title('x1(n)');

第 3 页

%======================8点======================= X1kk=myfft(x1n,8);

subplot 334;stem(k1/4,abs(X1kk),'.');

xlabel('\\omega/\\pi');ylabel('幅度');title('8点DFT[x1(n)]的幅频特性'); subplot 335;stem(k1/4,angle(X1kk),'.')

xlabel('\\omega/\\pi');ylabel('相位');title('8点DFT[x1(n)]的相频特性'); X1kk8=fft(x1n,8); subplot 336;mstem(X1kk8);

xlabel('\\omega/\\pi');ylabel('幅度');title('8点DFT理想幅频特性'); axis([0,2,0,4]);

%======================16点======================= X1k=myfft(x1n,16);

subplot 337;stem(k2/8,abs(X1k),'.');

xlabel('\\omega/\\pi');ylabel('幅度');title('16点DFT[x1(n)]的幅频特性'); subplot 338;stem(k2/8,angle(X1k),'.')

xlabel('\\omega/\\pi');ylabel('相位');title('16点DFT[x1(n)]的相频特性'); X1kk16=fft(x1n,16);subplot 339;mstem(X1kk16);

xlabel('\\omega/\\pi');ylabel('幅度');title('16点DFT理想幅频特性'); axis([0,2,0,4]);

运行结果如下：

第 4 页

第 5 页

第 6 页

6. 总结分析与结论

（1）对于时域连续信号x6(t)cos8tcos16tcos20t,取采样频率

Fs32Hz20Hz，采样点N=16，分辨率FFs/N2Hz，TpN/Fs0.5s为一个周期。

（2）自编DFT函数采用公式，对于有限长度序列x(n)，其长度为M，其N点离散傅里叶变换X(k)DFT[x(n)]x(n)Wn0N1knNk0,1,,N1，WNej2N，

N称为DFT的变换区间长度，NM。

（3）将DFT变换后幅度谱与相应的FFT变换所得幅度谱进行对比，谱线基本一致。而在模拟信号x6(t)的谱分析中，进行8点DFT与FFT变换时，由于采样频率为32Hz，x6(t)周期为0.5s，采样8点仅为半个周期的信号，因此得到的离散谱不能作为代表周期信号x6(t)的频谱。同理，x5(n)的8点DFT同样不能代表其频谱特性。

（4）由x1(n)的8点和16点DFT，分别是x1(n)的频谱函数的8点和16点采样；因为x3(n)x2((n3))3R8(n)，故x2(n)和x3(n)的8点DFT的模相等，但是

第 7 页

N=16时二者不满足循环移位关系，谱线不同；x4(n)周期为8，8点和16点DFT均得到正确的单一频率正弦波的频谱，即在0.25处有一根单一谱线，x5(n)周期为16，其16点DFT所得频谱为原频谱函数采样所得，在0.1250.50.625处有单一谱线；x6(t)cos8tcos16tcos20t，进行采样后DFT变换，可见在4、8、10Hz处有三根谱线，变换区间N=16,32时，观察时间Tp0.5s，1s，是x6(t)的整数周期，频谱正确，且N=32时频谱幅度是N=16时的2倍，可验证用DFT对中期序列谱分析的理论。

（5）误差产生原因：

①混叠现象：如果采样频率不满足Fs2fc，就会在fFs/2附近产生较大误差，发生频率混叠。

②栅栏效应：用FFT作频谱分析时，得到的是离散谱，而信号（周期信号除外）是连续谱，只有当N较大时，离散谱的包络才能逼近于连续谱，因此N要适当选择大一些。有限长序列尾部补零，无限长序列可增大截取长度及DFT变换区间长度。

③截断效应：使用窗函数(n)对无限长信号x(n)进行截取得到y(n)进行谱分析，会使原来的离散谱线向附近展宽，发生泄漏；主谱线两边形成很多旁瓣，引起谱间干扰。 ④频率分辨率不够。

（6）FFT做谱分析时的参数选择：

①周期信号的频谱是离散谱，用整数倍周期的长度作FFT，得到的离散谱才能代表周期信号的频谱。

②对模拟信号进行谱分析时，首先要按照采样定理将其变成时域离散信号。如果是模拟周期信号，应该选取整数倍周期的长度，经过采样后形成周期序列，按照周期序列的谱分析进行。各个参数满足：Fs2fc，谱分辨率

FFs/N1/TP1/nT。

③有限长序列可在尾部补零，无线长序列增大截取长度与DFT变换区间长度。加长原始信号截取长度/Tp可提高频率分辨率F。

第 8 页

实验二：用双线性变换法设计IIR数字滤波器 1.实验目的

（1）熟悉用双线性变换法设计IIR数字滤波器的原理与方法。 （2）掌握数字滤波器的计算机仿真方法。

2.实验环境

MATLAB R2016a。

3.实验报告要求

（1）简述实验原理及目的。

（2）由所绘图的频响特性曲线及设计过程，简述双线性变换法的特点。

4.实验原理

设计IIR数字滤波器一般采用间接法（脉冲响应不变法和双线性变换法），应用最广泛的是双线性变换法。其基本设计过程是：①先将给定的数字滤波器指标转换成过渡模拟滤波器的指标；②设计过渡模拟滤波器；③将过渡模拟滤波器的系统函数转换成数字滤波器的系统函数。MATLAB信号处理工具箱中的各种IIR数字滤波器设计函数都是采用双线性变换法。使用教材第6章介绍的滤波器设计函数butter、cheby1、cheby2和ellip进行IIR数字滤波器设计。

数字滤波器的实现调用MATLAB信号处理工具箱函数filter对给定的输入信号x(n)进行滤波得到滤波后输出信号y(n)。

5. 实验步骤

（1）复习关于butterworth模拟滤波器设计和用双线性法设计IIR数字滤波器的内容，按照例6.4.2，用双线性变换法设计数字滤波器系统函数H(z)。

例6.4.2中已求出满足本实验要求的数字滤波器系统函数：

0.0007378(1z1)6H(z)2.1(11.268z10.7051z2)(11.0106z10.3583z2)(10.9044z10.2155z2)Hk(z)k13式中，

A(12z1z2)Hk(z)k1,2,31Bkz1Ckz22.2

A=0.09036，B1=1.2686，B2=1.0106，B3=0.9044，C1=-0.7051，C2=-0.3583，C3=-0.2155。

由式2.1和2.2可以看出，滤波器H(z)由三个二阶滤波器H1(z)、H2(z)和

第 9 页

H3(z)级联组成。如图所示：

（2）编写滤波器仿真程序，计算H(z)对心电信号采样序列x(n)的响应序列y(n)。设yk(n)为第k级二阶滤波器Hk(z)的输出序列，yk1(n)为输入序列，由式2.2得到差分方程：

yk(n)Ayk1(n)2Ayk(n1)Ayk1(n2)Bkyk(n1)Ckyk(n2)2.3

当k=1时，yk1(n)x(n).所以H(z)对x(n)的总响应序列y(n)可以用顺序迭代算法得到。即以此对k=1,2,3，求解差分方程式2.3，最后得到y3(n)y(n)。仿真程序就是实现上述求解差分方程式和顺序迭代算法的通用程序。也可以直接调用Matlab的filter库函数实现仿真。

（3）运行仿真滤波程序，并绘图滤波器幅频响应特性曲线。

6. 实验内容

（1）用双线性变换法设计一个butterworth 低通IIR数字滤波器。设计指标参数为：在通带内频率低于0.2π时，最大衰减小于1dB，在阻带内[ 0.3π， π]的频率区间上，最小衰减大于15dB。

（2）打印出数字滤波器在频率区间 [0,π]上的幅频响应特性曲线。 ①确定数字低通的技术指标：

p0.2rad,p1dBs0.3rad，s15dB

②取T=1s，预畸变校正计算相应模拟低通的技术指标为：

p2tan2tan0.10.6498rad/s,p1dBT2

s2stan2tan0.151.0191rad/s,s15dBT2p③设计butterworth低通模拟滤波器，根据教材（6.2.18）式阶数N计算如下：

sps1.01911.568p0.649810s/101 ksp10.8751p/10101lgksp1g10.8751N5.3056lgsplg1.568

第 10 页

取N=6。将s和s代入（6.2.20）式，求得c=0.7662rad/s。这样保证阻带技术指标满足要求，通带指标有富余。

根据N=6，查表6.2.1得到归一化低通原型系统函数Ga(p)为

Ga(p)1 222(p0.5176p1)(p1.4142p1)(p1.9319p1)将ps/c代入Ga(p)，去归一化得到实际的Ha(s)为

Ha(s)0.2024

(s20.396s1)(s21.083s0.5871)(s21.480s0.5871)④用双线性变换法将Ha(s)转换成数字滤波器H(z)，即

0.0007378(1z1)6H(z)(11.268z10.7051z2)(11.0106z10.3583z2)(10.9044z10.2155z2)数字滤波器设计以及损耗函数绘制程序代码如下：

clc,clear; T=1;Fs=1/T; wpz=0.2;wsz=0.3;

wp=2*tan(wpz*pi/2);ws=2*tan(wsz*pi/2);rp=1;rs=15;%预畸变校正转换指标 [N,wc]=buttord(wp,ws,rp,rs,'s'); %设计过渡模拟滤波器，计算数字滤波器的阶数和3dB截止频率

[B,A]=butter(N,wc,'s'); %绘制模拟滤波器的损耗函数曲线 fk=0:1/512:1; wk=2*pi*fk; Hk=freqs(B,A,wk); subplot(1,2,1);

plot(fk,20*log10(abs(Hk))); grid on;

xlabel('f/Hz');ylabel('幅度(dB)');

axis([0,1,-100,5]);title('(a)模拟滤波器幅度特性'); [Bz,Az]=bilinear(B,A,Fs); %绘制数字滤波器的损耗函数曲线 Wk=0:pi/512:pi; subplot(1,2,2); Hk=freqz(Bz,Az,Wk); subplot(1,2,2);

plot(Wk/pi,20*log10(abs(Hk))); grid on;

第 11 页

xlabel('\\omega/\\pi');ylabel('幅度(dB)'); axis([0,1,-100,5]);title('(b)数字滤波器幅度特性');

（3）运用MATLAB产生两个正弦信号，信号频率为50Hz和400Hz，采样频率为1000Hz。两个正弦信号相叠加为输入信号y(t)。设计一滤波器，保留源信号中 50Hz的低频信号，对y(t)信号进行滤波。观察滤波前后信号的频谱特性，评价滤波器效果。

根据输入信号y(t)设计滤波器如下：

clc,clear;

wpd=0.2*pi;wsd=0.3*pi;%滤波器阻带与通带截止频率 Rp=1;As=15; %滤波器的通阻带衰减指标 Fs=1;T=1/Fs;

wp=(2/T)*tan(wpd/2); ws=(2/T)*tan(wsd/2);

[N,wc]=buttord(wp,ws,Rp,As,'s'); [z0,p0,k0]=buttap(N); ba=k0*poly(z0); aa=poly(p0);

[Ba1,Aa1]=lp2lp(ba,aa,wc);%变换为模拟低通滤波器 %用双线性变换法计算数字滤波器系数

[Bd,Ad]=bilinear(Ba1,Aa1,Fs) %双线性变换 %求数字系统的频率特性 [H,w]=freqz(Bd,Ad);

dbH=20*log10(abs(H)/max(abs(H))); %化为分贝值 subplot 321;plot(w/pi,abs(H));

grid on;ylabel('|H|');xlabel('\\omega/\\pi'); title('(a)数字滤波器幅度响应'); axis([0,1,0,1.1]);

subplot 322;plot(w/pi,dbH);

grid on; ylabel('dB');xlabel('\\omega/\\pi'); title('(b)数字滤波器幅度响应（dB)'); axis([0,1,-100,5]);

%====================================================== fs=1000; T=1/fs;Tp=0.1;t=0:T:Tp;%采样频率1000Hz x1=cos(2*pi*50*t); x2=sin(2*pi*400*t); y=x1+x2;

subplot 323;plot(t,y); title('(c)滤波前输入信号y(t)');

sf=filter(Bd,Ad,y);%用所设计的滤波器进行滤波 subplot 324;plot(t,sf);

title('(d)滤波后信号y1(t)'); axis([0,0.1,-1,1]); N=1024;T=1/fs;

第 12 页

k=0:N-1; fw=fft(y,N);

subplot 325; plot(k/512,abs(fw));

xlabel('\\omega/\\pi');title('(f)滤波前频谱'); fw1=fft(sf,N); subplot 326;

plot(k/512,abs(fw1));

xlabel('\\omega/\\pi');title('(g)滤波后频谱');

7. 总结分析与结论

（1）双线性变换法设计 butterworth 低通 IIR 数字滤波器的实验图象如图：

（2）根据输入信号y(t)cos(800t)cos(100t)设计低通滤波器幅度响应如图：

输入信号y(t)滤波前后及频率响应曲线如下：

第 13 页

滤波效果评价：

由滤波前后频谱可以看出，设计的低通滤波器通带截止频率wp=0.2π，通带最大衰减1dB,阻带截止频率为0.3π阻带最小衰减15dB。

滤波后400Hz的高频信号谱线消失，只留下50Hz的低频信号谱线，滤波效果较好；但是滤波后得到的信号y1(t)样本初始第一个周期内略有失真，且发生π/4的相位偏移。滤波总体效果符合要求。

双线性变换法的特点：

第 14 页

①模拟滤波器经过双线性变换后，不存在频率特性的混叠失真，因而对模拟滤波的频率响应函数Ha(s)无限带要求，而且能够直接用于设计低通、高通、带通、带阻等各种类型的数字滤波器。

②模拟频率ω与数字频率Ω之间的非线性关系，使数字滤波器响应曲线不能保真地模仿模拟滤波器的频响曲线形状，在高频处有较大的相频特性失真。事先进行频率预畸变，这种非线性关系不会使所设计的数字滤波器的幅频特性受到影响。

③双线性变换法无需将模拟系统函数进行部分因式分解。

第 15 页

实验三：语音信号谱分析及去噪处理

1.实验目的

（1）通过对实际采集的语音信号进行分析和处理，获得数字信号处理实际应用的感性认识。

（2）掌握数字信号谱分析的知识。

（3）掌握数字滤波器设计的知识，并通过对语音信号的去噪处理，获得数字滤波器实际应用的感性知识。

2.实验环境

MATLAB R2016a。

3.实验报告要求

（1）简述实验原理及目的。

（2）简述数字滤波器指标设定的计算。

（3）对比去噪前后的语音信号，讨论数字滤波器的滤波作用。

4.实验内容

（1）用麦克风自行采集两段语音信号[高频噪声、人声+高频噪声]（.wav 格式）。 （2）通过Matlab读入 采集信号，观察其采样频率，并绘图采样信号。 （4）通过Matlab对语音信号进行谱分析，分析出噪声的频带。

（5）设计一滤波器，对叠加入噪声的语音信号进行去噪处理。绘图并发声去噪后的信号。

5.实验步骤

（1）利用麦克风采集一段5s以内的语音信号。利用格式工厂软件对语音信号进行预处理。通常语音信号为单声道，采样频率为8000Hz，语音信号为.wav格式。 （2）通过Matlab读入语音信号及其采样频率（使用Matlab库函数wavread），在Matlab软件的workspace工作平台上观察读入的语音信号，在Matlab中，对入的语音信号为一维矩阵。应注意，库函数wavread自动将语音信号幅度归一化[-1,1]区间范围。使用Matlab库函数plot绘图语音信号，并使用库函数sound 发音语音信号。

（3）分析噪声的频谱。在这里进行谱分析的目的，是了解噪声信号的频谱特性，为去噪滤波器的技术指标提供依据。叠加入噪声的语音信号，如下图所示：

第 16 页

（4）通过Matlab对语音信号进行谱分析。应注意，对信号进行谱分析，在实验一中已经详细介绍过。在这里进行谱分析的目的，是了解本段语音信号的频谱特性，为去噪滤波器的技术指标提供依据。

（5）根据语音信号及噪声信号的频谱特性，自行设计一滤波器，对叠加入噪声的语音信号进行去噪处理。最后绘图并发声去噪后的信号。应注意，数字滤波器的实际应考虑实际需求，合理制定滤波器的技术指标。

实验代码：

clc,clear;

%===================语音信号=======================

[x,Fs]=audioread('Fire.wav');%调用音频文件，采样值放在x中，fs为采样频率 x=x(:,1); %取矩阵x的第一列赋值到x矩阵中 FS=length(x); % 矩阵x的长度 f=0:Fs/FS:(FS-1)*Fs/FS; %生成一个一维数组赋给f X=fft(x,65536); %对信号做FFT变换 2^16=65536 magX=abs(X); %取X的幅值给magX t=(0:FS-1)/Fs;

subplot 221;plot(t,x); %绘制原始语音信号的时域波形图 title('原始语音信号时域波形图'); axis([0,5,-1,1]);

xlabel('时间/n'); ylabel('幅值/n');grid on; f=Fs*(0:511)/1024; %生成一个一维数组赋给f subplot 223; plot(0:2*pi/65535:2*pi,magX); title('原始语音信号频谱图');grid on; %sound(x,fs)

%====================高频噪声=========================== nt=2*rand(1,FS)-1;%随机噪声nt fp=2500;fs=3000;Rp=0.1;As=70;

fb=[fp,fs];m=[0,1]; %计算remezord函数所需参数f,m,dev dev=[10^(-As/20),(10^(Rp/20)-1)/(10^(Rp/20)+1)]; [n,fo,mo,W]=remezord(fb,m,dev,Fs);%确定remez函数所需参数

hn=remez(n,fo,mo,W); %调用remez函数进行设计，用于滤除噪声nt中的低频成分 yt=2*filter(hn,1,10*nt)/25; %滤除随机噪声中低频成分，生成高通噪声yt Y=fft(yt,65536); %对信号做FFT变换 2^16=65536 magY=2*abs(Y); %取X的幅值给magX

subplot 222;plot(t,yt); %绘制原始语音信号的时域波形图 title('高频噪声信号时域波形图'); axis([0,5,-1,1]);

xlabel('时间/n'); ylabel('幅值/n');grid on; f=fs*(0:511)/1024; %生成一个一维数组赋给f

第 17 页

subplot 224; plot(magY); title('高频噪声信号频谱图'); grid on;

%======================叠加噪声=============================== y=x'+yt;Y=fft(y,65536); %对信号做变换 2^16=65536 magY=abs(Y); %取Y的幅值给magY t=(0:FS-1)/Fs;

figure(2);subplot 241;plot(t,y); %绘制原始语音信号的时域波形图 title('滤波前信号时域波形图'); axis([0,5,-1,1]);

xlabel('时间/n'); ylabel('幅值/n');grid on; f=Fs*(0:511)/1024; %生成一个一维数组赋给f subplot 245; plot(0:2*pi/65535:2*pi,magY); title('滤波前信号频谱图'); grid on; %sound(x,fs)

%==============设计低通滤波器fp=2500，fs=2700============= N=40000;fp=2500;fs=2700;Rp=0.2;As=60; %(1)窗函数法设计滤波器

wc=(fp+fs)/Fs; %理想低通滤波器截止频率（关于π归一化） B=2*pi*(fs-fp)/Fs;%过渡带宽度指标 Nb=ceil(11*pi/B); %blackman窗长度N hn=fir1(Nb-1,wc,blackman(Nb));

Hw=abs(fft(hn,1024));%求设计的滤波器频率特性 figure(3);subplot 311;

myplot(hn,1);title('窗函数法滤波器损耗函数曲线');

y1=fftfilt(hn,y,N);Y1=fft(y1,65536); %对信号做FFT变换 2^16=65536 magY1=abs(Y1); %取Y1的幅值给magY1

figure(2)subplot 242;plot(t,y1); %绘制原始语音信号的时域波形图 title('滤波后语音信号时域波形图'); axis([0,5,-1,1]);

xlabel('时间/n'); ylabel('幅值/n');grid on; f=Fs*(0:511)/1024; %生成一个一维数组赋给f subplot 246; plot(0:2*pi/65535:2*pi,magY1); title('滤波后语音信号频谱图');grid on; %sound(x,fs) %等波纹最佳逼近法设计滤波器

fb=[fp,fs];m=[1,0];%确定remez函数所需参数

dev=[(10^(Rp/20)-1)/(10^(Rp/20)+1),10^(-As/20)]; [Ne,fo,mo,W]=remezord(fb,m,dev,Fs);%确定remez函数所需参数 hn=remez(Ne,fo,mo,W); %调用remez函数进行设计 figure(3);subplot 312;

myplot(hn,1);title('等波纹最佳逼近法滤波器损耗函数曲线'); Hw=abs(fft(hn,1024)); %求设计的滤波器频率特性 y2=fftfilt(hn,x,N); %调用函数fftfilt对x滤波 Y2=fft(y2,65536); %对信号做FFT变换 2^16=65536 magY2=abs(Y2); %取Y2的幅值给magY2

figure(2)subplot 243;plot(t,y2); %绘制原始语音信号的时域波形图 title('滤波后语音信号时域波形图'); axis([0,5,-1,1]);

xlabel('时间/n'); ylabel('幅值/n');grid on;

第 18 页

subplot 247; plot(0:2*pi/65535:2*pi,magY2); title('滤波后语音信号频谱图');grid on; %(3)双线性变换法设计滤波器。

wpd=fp/Fs*2*pi;wsd=fs/Fs*2*pi;%滤波器阻带与通带截止频率 FS=1;T=1/FS;

wp=(2/T)*tan(wpd/2);ws=(2/T)*tan(wsd/2); [N,wc]=buttord(wp,ws,Rp,As,'s'); [z0,p0,k0]=buttap(N);

ba=k0*poly(z0); aa=poly(p0); [Ba1,Aa1]=lp2lp(ba,aa,wc);%变换为模拟低通滤波器 [Bd,Ad]=bilinear(Ba1,Aa1,FS); %双线性变换 [H,w]=freqz(Bd,Ad); figure(3);subplot 313;

myplot(Bd,Ad);title('双线性变换法滤波器损耗函数曲线');

y3=filter(Bd,Ad,x); Y3=fft(y3,65536); %对信号做FFT变换 2^16=65536 magY3=abs(Y3); %取Y3的幅值给magY3

figure(2);subplot 244;plot(t,y3); %绘制原始语音信号的时域波形图 title('滤波后语音信号时域波形图'); axis([0,5,-1,1]);

xlabel('时间/n'); ylabel('幅值/n');grid on; subplot 248; plot(0:2*pi/65535:2*pi,magY3); title('滤波后语音信号频谱图');grid on;

6.实验结果分析与结论

（1）因为MATLAB R2016A版本没有wavread函数，故使用audioread('Fire.wav')读入语音信号。

语音信号读入后赋值给x，由其频谱特性可以看出，所分析的语音信号频谱

分布并不单一，有许多不同频率的成分，但是在大于2500Hz处谱线的幅度很小，即其在原语音信号中所占的成分并不多。

所以在设计产生高频噪声信号时，先产生所有频率的随机噪声信号nt，再使用高通数字滤波器进行滤波产生高频噪声信号yt。由上面的分析过程，所设

第 19 页

计高通滤波器的指标为fp=2500;fs=3000;Rp=0.1;As=70。

叠加处理前的语音与高频噪声信号及各自频谱如下：

直接使语音信号x与高频噪声信号yt相叠加产生滤波处理前信号y，如图。 （2）使用三种方法设计低通滤波器：

①窗函数法设计FIR低通滤波器；

②等波纹最佳逼近法设计FIR低通滤波器； ③双线性法设计巴特沃斯低通滤波器。

其中滤波器指标开始时设置为fp=2500;fs=3000;Rp=0.2;As=60，具体计算过程在程序中已经列出。发现设计出来的滤波器滤波效果并不很好，频谱图里噪声成分残留较多；后调整为fp=2500;fs=2700;Rp=0.2;As=60;发现如此三个数字滤波器的滤波性能都大大提高。

数字滤波器损耗函数曲线如下：

第 20 页

对叠加后的语音信号分别进行滤波处理，得到滤波后信号时域波形和频谱特性如下：

与原始语音信号时域波形和频谱图线相对比，三种滤波器滤波效果都很好。

（3）使用sound(y,8000)函数对语音信号x、y、y1、y2、y3分别进行发声，发现三种数字滤波器滤波后的语音信号y1、y2、y3都出现一定程度的损耗，出现一定程度的失真。

其中，使用窗函数法设计的滤波器声音信号中还存有少量噪声，并没有完全滤除高频噪声；双线性变换法设计的巴特沃斯低通滤波器滤波后信号y3也存有一些噪声，但是比窗函数法滤波器的效果要好；等波纹逼近

法设计的滤波器滤波后信号y2基本听不出噪声的存在，其滤波效果是三个低通滤波器中最好的。

（4）综上分析，所设计的三个数字滤波器都达到了设计要求，能够有效地从噪声中提取信号。

第 21 页

实验心得

参考文献

[1] 数字信号处理[M] 丁玉美 高西全等 西安电子科技大学出版社.2008年 [2] 数字信号处理学习指导[M] 丁玉美 高西全等 西安电子科技大学出版社.2009年

第 22 页