数字信号处理课程设计报告语音信号的数字滤波处理大学论文.docx
《数字信号处理课程设计报告语音信号的数字滤波处理大学论文.docx》由会员分享，可在线阅读，更多相关《数字信号处理课程设计报告语音信号的数字滤波处理大学论文.docx（31页珍藏版）》请在冰豆网上搜索。
数字信号处理课程设计报告语音信号的数字滤波处理大学论文
信号处理综合设计
设计说明书
设计项目：
语音信号的数字滤波处理
项目完成人：
姓名：
学号：
专业班级：
电子信息工程二班
指导教师：
提交日期：
2017.1.7
机电工程学院电子信息工程
1.设计基本原理
1.1课题研究的背景
1）熟悉离散信号和系统的时域特性。
2）掌握数字信号处理的基本概念，基本理论和基本方法。
3）掌握序列快速傅里叶变换方法。
4）练习对MATLAB的使用，掌握MATLAB的程序设计方法。
5）学习利用MATLAB对语音信号进行频谱分析。
6）学习设计巴特沃斯滤波器以及采用hamming窗设计FIR数字低通，高通和带通滤波器的方法。
1.2课题研究意义
通过本次课程设计，综合运用数字信号处理技术课程和其他有关先修课程的理论和生产实际知识去分析和解决具体问题，并使所学知识得到进一步巩固、深化和发展。
初步培养学生对工程设计的独立工作能力，学习设计的一般方法。
通过课程设计树立正确的设计思想，提高学生分析问题、解决问题的能力。
进行设计基本技能的训练，如查阅设计资料和手册、程序的设计、调试，了解语音信号的产生、采集、能绘制语音信号的频率响应曲线及频谱图，学会用MATLAB对语音信号进行分析和处理，学习设计滤波器的一般方法，通过课程设计树立正确的设计思想，提高分析问题、解决问题的能力，掌握用滤波器去除语音信号噪声的方法，观察去噪前后的语音信号，进行设计基本技能的训练，如查阅设计资料和手册、程序的设计、调试等。
设计任务
2.1课题设计的任务
（1）熟练掌握数字信号处理的基本概念基本理论。
（2）熟悉离散信号和离散系统的时域特性。
（3）掌握序列快速傅里叶变换（FFT）方法。
（4）掌握滤波器的几种网络结构。
（5）学会使用软件MATLAB，掌握其程序设计方法。
（6）掌握利用MATLAB对语音信号进行频谱分析的方法。
2.2课题研究的内容
预习题部分
1、设计卷积运算的演示程序:
1可输入任意两个序列x1（n）、x2（n），指定x1（n）为自己的学号，例如x1（n）={2,0,1,0,5,7,0,5,0,2,0,1}。
x2（n）的内容和长度自选。
例如x2（n）={0.5,2.43,6.17,12.93,22.17,32.25,40.88,45.87,45.87,40.88,32.25,22.17,12.93,6.17,2.43,1.007}。
2分别动态演示两个序列进行线性卷积x1（n）*x2（n）和圆周卷积x1（n）⊙x2（n）的过程；要求分别动态演示翻转、移位、乘积、求和的过程；
3圆周卷积默认使用两个序列中的最大长度，但卷积前可以指定卷积长度N用以进行混叠分析；
4改变圆周卷积长度N，根据实验结果分析两类卷积的关系；
5在计算机操作系统中选一段声音文件（XP系统在“C:
\WINDOWS\Media”），读取文件取10ms的声音数据产生时域序列x1（n），序列内容自定义。
利用x2（n）={0.5,2.43,6.17,12.93,22.17,32.25,40.88,45.87,45.87,40.88,32.25,22.17,12.93,6.17,2.43,1.007}。
利用FFT实现快速卷积，验证时域卷积定理，并与直接卷积进行效率对比（验证时采用Matlab子函数）。
2、编写程序演示采样定理（时域采样、频谱周期延拓），同时演示采样频率小于2fc时，产生的混叠效应：
1对下面连续信号进行采样：
学号，A为幅度因子，a为衰减因子，
为模拟角频率，其中n为学号（例如，贺娜同学n=201）；
2要求输入采样频率fs（根据程序处理需要指定范围）后，在时域演示信号波形、采样脉冲及采样后信号；在频域演示不同采样频率下对应信号的频谱。
设计题部分
数字滤波器设计及其应用
1利用Windows下的录音机或其他软件，进行语音信号的采集（*.wav）；
2语音信号的频谱分析，画出采样后语音信号的时域波形和频谱图；
3产生噪声信号并加到语音信号中，得到被污染的语音信号，并回放语音信号；
4污染信号的频谱分析，画出被污染的语音信号时域波形和频谱；
5根据有关的频谱特性，采用间接法设计IIR数字滤波器，并画出相应滤波器的幅频、相频图（设计3个IIR滤波器）
a．模拟滤波器类型：
巴特沃斯滤波器（低通、带通、高通）
b．总体要求：
Matlab原程序+仿真波形+技术指标
6根据有关的频谱特性，采用直接法设计FIR数字滤波器，并画出相应滤波器的幅频、相频图（设计3个FIR滤波器）
a．滤波器类型：
汉宁（Hanning）窗（低通、带通、高通）
b．总体要求：
Matlab原程序+仿真波形+技术指标+窗函数
7用自己设计的这些滤波器分别对被不同噪声污染的信号进行滤波；
8分析得到信号的频谱，画出滤波后信号的时域波形和频谱，并对滤波前后的信号进行对比，分析信号的变化；
回放语音信号。
设计提示
1．语音信号的采集
要求利用windows下的录音机（开始—程序—附件—娱乐—录音机，文件—属性—立即转换—8000Hz，8位，单声道）录制一段自己的话音，或者采用Windows自带的声音文件（默认为22050Hz），时间控制在几秒左右。
然后在MATLAB软件平台下，利用函数wavread对语音信号进行采样，记住采样频率和采样点数。
通过wavread函数的使用，要求理解采样频率、采样位数等概念。
wavread函数调用格式：
y=wavread（file），读取file所规定的wav文件，返回采样值放在向量y中。
[y,fs,nbits]=wavread（file），采样值放在向量y中，fs表示采样频率（Hz），nbits表示采样位数。
y=wavread（file,N），读取前N点的采样值放在向量y中。
y=wavread（file,[N1,N2]），读取从N1点到N2点的采样值放在向量y中。
2．语音信号的频谱分析
要求首先画出语音信号的时域波形；然后对语音信号进行频谱分析，在MATLAB中，可以利用函数fft对信号进行快速付立叶变换，得到信号的频谱特性；从而加深对频谱特性的理解。
3．设计数字滤波器和画出频率响应
根据语音信号的特点给出有关滤波器的性能指标，例如：
1）低通滤波器性能指标，fp=1000Hz，fc=1200Hz，As=100dB，Ap=1dB；2）高通滤波器性能指标，fc=2800Hz，fp=3000Hz，As=100dB，Ap=1dB；3）带通滤波器性能指标，fp1=1200Hz，fp2=3000Hz，fc1=1000Hz，fc2=3200Hz，As=100dB，Ap=1dB。
4.回放语音信号
在Matlab中,函数sound可以对声音进行回放。
其调用格式：
sound（x,fs,bits）;可以感觉滤波前后的声音有变化。
5.语音的反转
将女生的一段声音反转成男生，或者将男生的声音反转成女生。
2.3课题研究的指标
在设计一个滤波器之前，必须首先根据工程实际的需要确定滤波器的技术指标。
在很多实际应用中，数字滤波器常常被用来实现选频操作。
因此，指标的形式一般在频域中给出幅度和相位响应。
幅度指标主要以两种方式给出。
第一种是绝对指标。
它提供对幅度响应函数的要求，一般应用于FIR滤波器的设计。
第二种指标是相对指标。
它以分贝值的形式给出要求。
在工程实际中，这种指标最受欢迎。
对于相位响应指标形式，通常希望系统在通频带中人有线性相位。
运用线性相位响应指标进行滤波器设计具有如下优点：
①只包含实数算法，不涉及复数运算；②不存在延迟失真，只有固定数量的延迟；③长度为N的滤波器（阶数为N-1），计算量为N/2数量级。
因此，本文中滤波器的设计就以线性相位FIR滤波器的设计为例。
预习部分练习
1.线性卷积
n1=1:
40;
n2=1:
40;
x1=[zeros（1,12）,2,0,1,4,2,2,1,0,1,2,0,0,zeros（1,16）];x2=[zeros（1,12）,2,4,5,6,2,7,4,8,9,6,4,8,3,5,7,8,zeros（1,12）];
x11=zeros（size（x1））;
subplot（3,1,1）
stem（n1,x1）
subplot（3,1,2）
stem（n2,x2）
pause
（1）;
fori=1:
24
x11（i）=x1（25-i）;
end
fori=25:
40
x11（i）=0;
end
subplot（3,1,1）
stem（n1,x11）;
pause
（1）
x3=zeros（size（x11））;
fori=1:
28
x11=[zeros（1,1）,x11（1:
（40-1））];
x3（i）=sum（x11.*x2）;
subplot（3,1,1）
stem（n1,x11）
subplot（3,1,3）
x4=[zeros（1,12）,x3（1:
i）,zeros（1,28-i）];
stem（n1,x4）
pause
（1）
end
圆周卷积
x=[2,0,1,4,2,2,1,0,1,2,0,0];
y=[5,4,3,6,1,7,4,8,4,5,4,8,3,8,7,8];
subplot（2,1,1）
stem（[x,zeros（1,28）]）;
N1=length（x）;
N2=length（y）;
N=N1;
ifN2>N1
N=N2;
end
M=input（'请输入圆周卷积长度\n'）;
ifM>=N
N=M;
end
ifN==N1
x=x;
else
x=[x,zeros（1,N-N1）];
end
X=zeros（N,N）;
x1=zeros（1,N）;
pause
（2）
forn=1:
N
ifn==1
x1（n）=x（n）;
else
x1（n）=x（N-n+2）;
end
end
X（1,1:
N）=x1;
fori=2:
N
pause
（1）
forn=1:
N
ifn==1
X（i,n）=X（i-1,N）;
else
X（i,n）=X（i-1,n-1）;
end
end
subplot（2,1,1）
stem（X（i,1:
N））;
end
ifN2y=[y,zeros（1,N-N2）];
end
Y=y';
W=X*Y;
W=W';
n=1:
N;
subplot（2,1,2）
stem（n,W）
分析总结：
由以上分析可知，当循环卷积区间长度大于28时，不会出现混叠显现；当循环卷积长度小于28时，会出现混叠现象从而失真。
利用FFT实现快速卷积
x1=wavread（'Alarm01.wav',2000）;
x1=x1';x1=x1（1,:
）;
x2=[12.436.1712.9322.1732.2540.8845.8745.8740.8832.2522.1712.936.172.431];
subplot（7,1,1）;
stem（x1）;
ylabel（'x1（n）'）;
title（'x1（n）'）;
subplot（7,1,2）;
stem（x2）;
ylabel（'x2（n）'）;
title（'x2（n）'）;
y=conv（x1,x2）;
subplot（7,1,3）;
stem（y）;
ylabel（'y'）;
title（'x1（n）与x2（n）的卷积'）;
N1=length（x1）;
N2=length（x2）;
N=N1+N2-1;
X1=fft（x1,N）;X2=fft（x2,N）;
subplot（7,1,4）;
stem（X1）;
ylabel（'X1'）;
title（'x1（n）的N点DFT'）;
subplot（7,1,5）;
stem（X2）;
ylabel（'X2'）;
title（'x2（n）的N点DFT'）;
Y1=X1.*X2;
subplot（7,1,6）;
stem（Y1）;
ylabel（'Y1'）;
title（'X1与X2相乘的结果'）;
Y2=ifft（Y1）;
subplot（7,1,7）;
stem（Y2）;
ylabel（'Y2'）;
title（'Y1的IDFT结果'）;
fs=input（'请输入采样频率:
'）;
n=input（'请输入4位学号:
'）;
t=0:
1/fs:
1;
t1=0:
0.0005:
1;
A=n;
a=sqrt
（2）*n*pi;
w=a;
x=A*exp（-a*t）.*sin（w*t）;
x1=A*exp（-a*t1）.*sin（w*t1）;
subplot（5,1,1）
plot（x1）;
xn=ones（1,length（t））;
subplot（512）
stem（xn）
yn=xn.*x;
subplot（5,1,3）
stem（yn）;
y=fft（yn,fs）;
f=fs*（0:
fs/2）/fs;
subplot（5,1,4）
plot（f,y（1:
fs/2+1））;
yn=ifft（y）;
subplot（5,1,5）
plot（yn）
编写程序演示采样定理（时域采样、频谱周期延拓），同时演示采样频率小于2fc时，产生的混叠效应（n=14）：
t=0:
0.0005:
0.3;
xa=14*exp（-14*sqrt
（2）*pi*t）.*sin（14*sqrt
（2）*pi*t）;
subplot（6,1,1）;plot（t,xa）;title（'originxa'）
t=0.015;
n=0:
t:
0.3;
pt=ones（0.3/t+1）;
subplot（6,1,2）;stem（n,pt）;title（'采样脉冲'）;
xaa=14*exp（-14*sqrt
（2）*pi*n）.*sin（14*sqrt
（2）*pi*n）;
subplot（6,1,3）;stem（n,xaa）;title（'采样后信号'）;
y1=fft（xa,4096）;
subplot（6,1,4）;plot（abs（y1））;title（'采样频率2kHz'）;
y2=fft（xaa,4096）;
subplot（6,1,5）;plot（abs（y2））;title（'采样频率66Hz'）;
t=0.005;
n=0:
t:
0.3;
xa1=14*exp（-14*sqrt
（2）*pi*n）.*sin（14*sqrt
（2）*pi*n）;
y2=fft（xa1,4096）;
subplot（6,1,6）;plot（abs（y2））;title（'采样频率200Hz'）;
结果分析总结：
当序列的长度为2^N时，可用FFT变换，再用两者频域相乘，逆变换成时域，计算速度比直接计算卷积快的多。
当循环区间长度L>=N+M-1时，线性卷积与循环卷积相等，且无时域混叠现象。
滤波器的设计
除了低通数字滤波器之外，实际中还常常需要高通、带通及带阻数字滤波器。
这三种数字滤波器的设计步骤如下：
（1）将数字滤波器H（z）的技术指标ωp和ωs，通过Ω=tan（ω/2）转变为模拟滤波器G（s）的技术指标Ωp和Ωs，作归一化处理后，得到ηp=1，ηs=Ωs/Ωp；
（2）化解为模拟原型滤波器G（s）的技术指标；
（3）设计模拟原型滤波器G（p）；
（4）将G（p）转换为模拟滤波器的转移函数G（s）；
（5）将G（s）转换成数字滤波器的转移函数H（z）s=（z–1）（z+1）。
所谓原型滤波器是指归一化的低通滤波器。
本节主要讨论通过IIR数字滤波器的原型转换设计法和IIR数字滤波器的直接设计方法来设计数字高通、带通及带阻滤波器其转换方法主要有3种一是直接由模拟低通滤波器转换成数字高通、带通或带阻滤波器二是先由模拟低通滤波器转换成模拟高通、带通或带阻滤波器然后再把它转换成相应的数字滤波器三是将模拟低通滤波器先转换成数字低通滤波器再通过变量代换变换成高通、带通或带阻滤波器。
语音信号的频谱范围主要为22050HZ左右，并且在50000HZ左右有一个小信号，因此，在设计低通滤波器时，应把噪声频谱设定在50000HZ以上，这样，通过低通滤波器，即可滤除噪声信号从而还原语音信号；在设计高通滤波器时，应把噪声设定在22050HZ以内，以通过高通滤波器滤除低频的噪声信号，从而还原相对频率较高的语音信号；在设计带通滤波器时，可把噪声设计在低于22050HZ或高于50000HZ的频谱上，已通过带通滤波器还原带通范围内的语音信号。
巴特沃思滤波器低通滤波器
%读入音频信号
[x,FS,bits]=wavread（'Alarm01.wav'）;
x=x（:
1）;
figure
（1）;
subplot（2,1,1）;
plot（x）;
sound（x,FS,bits）;
title（'语音信号时域波形图'）
y=fft（x,3260）;
f=（FS/1630）*[1:
1630];
subplot（2,1,2）;
plot（f（1:
1630）,abs（y（1:
1630）））;
title（'语音信号频谱图'）;
%产生噪声
t=0:
length（x）-1;
zs=0.05*cos（2*pi*6000*t/22050）;
zs0=0.05*cos（2*pi*6000*t/22050000）;
figure
（2）;
subplot（2,1,1）
plot（zs0）
title（'噪声信号波形'）;
zs1=fft（zs,1200）;
%sound（zs,FS,bits）;
subplot（2,1,2）
plot（f（1:
600）,abs（zs1（1:
600）））;
title（'噪声信号频谱'）;
%tjzs
x1=x+zs';
sound（x1,FS,bits）;%回放加入噪声后的语音
y1=fft（x1,1200）;%用傅里叶转为频域
figure（3）;
subplot（2,1,1）;plot（x1）;
title（'加入噪声后的信号波形'）;
subplot（2,1,2）;
plot（f（1:
600）,abs（y1（1:
600）））;
title（'加入噪声后的信号频谱'）;
%巴特沃斯低通滤波器
fp=3000;fs=3500;Fs=22050;
rp=1;rs=10;
wp=2*pi*fp/Fs;
ws=2*pi*fs/Fs;
Fs1=1;
wap=2*tan（wp/2）;
was=2*tan（ws/2）;
[N,wc]=buttord（wap,was,rp,rs,'s'）;
[B,A]=butter（N,wc,'s'）;
[Bz,Az]=bilinear（B,A,Fs1）;
figure（4）;
[h,w]=freqz（Bz,Az,512,Fs1*22050）;
plot（w,abs（h））;
title（'巴特沃斯低通滤波器'）;
xlabel（'频率（HZ）'）;ylabel（'耗损（dB）'）;
gridon;
yd=filter（Bz,Az,x1）;
figure（5）;
subplot（2,1,1）;plot（yd）;title（'滤波后信号波形'）;
ydd=fft（yd,1200）;
subplot（2,1,2）;plot（f（1:
600）,abs（ydd（1:
600）））;
title（'滤波后信号频谱'）;
sound（yd,FS,bits）
巴特沃斯带通滤波器
%语音信号
[x,fs,bits]=wavread（'Alarm01.wav'）;%输入参数为文件的全路径和文件名，输出的第一个参数
%sound（x,fs,bits）;%按指定的采样率和每样本编码位数回放
N=length（x）;%计算信号x的长度
fn=5000;%单频噪声频率
t=0:
1/fs:
（N-1）/fs;%计算时间范围，样本数除以采样频率
x=x（:
1）';%将双声道转为单声道
%加入噪声
y=x+0.01*sin（fn*2*pi*t）;%加噪声
sound（y,fs,bits）;%应该可以明显听出有尖锐的单频啸叫声
X=abs（fft（x））;Y=abs（fft（y））;%对原始信号和加噪信号进行fft变换,取幅度谱
X=X（1:
N/2）;Y=Y（1:
N/2）;%截取前半部分
deltaf=fs/N;%计算频谱的谱线间隔
f=0:
deltaf:
fs/2-deltaf;%计算频谱频率范围
figure
（1）;%新建图形
subplot（2,2,1）;plot（x）;%原始语音信号幅度图
xlabel（'timen'）;%定义X坐标
ylabel（'幅度n'）;%定义Y坐标
title（'原始语音信号shiyishi'）;%定义主题
axis（[0120000-1.11.1]）;gridon%设定X轴和Y轴范围
subplot（2,2,2）;%定义第二个子图
plot（f,X）;%原始语音信号幅度谱图
title（'语音信号幅度谱图shiyishi'）;%定义主题
xlabel（'频率（单位：
HZ）'）;%定义Y坐标
ylabel（'幅度谱'）;%定义X坐标
axis（[080000600]）;gridon%设定X轴和Y轴范围
subplot（2,2,3）;%定义第三个子图
plot（y）;%加噪声后的语音信号图
axis（[0120000-1.00.9]）;gridon%设定X轴和Y轴范围
xlabel（'timen'）;
ylabel（'幅度n'）;
title（'加入干扰后的语音信号'）;
subplot（'2,2,4'）;
plot（f,Y）;%加噪声后的语音信号幅度谱图>>axis（[080000700]）;gridon
xlabel（'频率（单位：
HZ）'）;ylabel（'幅度谱'）;
title（'加入干扰后信号幅度谱图'）;
%巴特沃斯带通滤波器
fsl=300;fpl=800;fpu=1200;fsu=2200;fs_1=22050;
rp=1;rs=10;
wp=2*pi*[fpl,fpu]/fs;
ws=2*pi*[fsl,fsu]/fs;
Fs=fs/fs;
wap=2*tan（wp/2）;was=2*tan（ws/2）;
[N,wc]=buttord（wap,was,rp,rs,'s'）;
[B,A]=butter（N,wc,'s'）;
[Bz,Az]=bilinear（B,A,Fs）;
[h,w]=freqz（Bz,Az,512,Fs*44100）;
figure
（2）;
subplot（311）;
plot（w/pi,abs（h））;
title（'巴特沃斯带通滤波器-频率响应'）;
xlabel（'频率（HZ）'）;ylabel（'耗损（dB）'）;
gridon;
yd=filter（Bz,Az,y）;
subplot（312）;
升级会员 