MATLAB实现数字信号处理.docx
《MATLAB实现数字信号处理.docx》由会员分享，可在线阅读，更多相关《MATLAB实现数字信号处理.docx（20页珍藏版）》请在冰豆网上搜索。
MATLAB实现数字信号处理
《数字信号处理》课程设计实例：
声音信号的处理
一．摘要：
这次课程设计的主要目的是综合运用本课程的理论知识进行频谱分析以及滤波器设计，通过理论推导得出相应结论，并利用MATLAB或者ＤＳＰ开发系统作为工具进行实现，从而复习巩固课堂所学的理论知识，提高对所学知识的综合应用能力，并从实践上初步实现对数字信号的处理。
通过对声音的采样，将声音采样后的频谱与滤波。
MATLAB全称是MatrixLaboratory，是一种功能强大、效率高、交互性好的数值和可视化计算机高级语言，它将数值分析、矩阵运算、信号处理和图形显示有机地融合为一体，形成了一个极其方便、用户界面友好的操作环境。
。
经过多年的发展，已经发展成为一种功能全面的软件，几乎可以解决科学计算中所有问题。
MATLAB软件还提供了非常广泛和灵活的用于处理数据集的数组运算功能。
在本次课程设计中，主要通过MATLAB来编程对语音信号处理与滤波，设计滤波器来处理数字信号并对其进行分析。
二．课程设计目的：
综合运用本课程的理论知识进行频谱分析以及滤波器设计，通过理论推导得出相应结论，并利用MATLAB作为工具进行实现，从而复习巩固课堂所学的理论知识，提高对所学知识的综合应用能力，并从实践上初步实现对数字信号的处理。
三．设计内容：
内容：
录制一段个人自己的语音信号，并对录制的信号进行采样；画出采样后语音信号的时域波形和频谱图；给定滤波器的性能指标，采用窗函数法和双线性变换法设计滤波器，并画出滤波器的频率响应；然后用自己设计的滤波器对采集的信号进行滤波，画出滤波后信号的时域波形和频谱，并对滤波前后的信号进行对比，分析信号的变化；回放语音信号；换一个与你性别相异的人录制同样一段语音内容，分析两段内容相同的语音信号频谱之间有什么特点；再录制一段同样长时间的背景噪声叠加到你的语音信号中，分析叠加前后信号频谱的变化，设计一个合适的滤波器，能够把该噪声滤除。
四．设计原理：
4.1．语音信号的采集
熟悉并掌握MATLAB中有关声音（wave）录制、播放、存储和读取的函数，在MATLAB环境中，有关声音的函数有：
a：
y=wavrecord（N,fs,Dtype）;利用系统音频输入设备录音，以fs为采样频率，默认值为11025，即以11025HZ进行采样。
Dtype为采样数据的存储格式，用字符串指定，可以是：
‘double’、‘single’、’int16’、‘int8’其中只有int8是采用8位精度进行采样，其它三种都是16位采样结果转换为指定的MATLAB数据；
b：
wavplay（y,fs）；利用系统音频输出设备播放，以fs为播放频率，播放语音信号y；
c：
wavwrite（（y,fs,wavfile）;创建音频文件；
d：
y=wavread（file）;读取音频文件；
关于声音的函数还有sound（）；soundsc（）；等。
4.2滤波器：
4．21.IIR滤波器原理
冲激响应不变法是使数字滤波器在时域上模拟滤波器，但是它们的缺点是产生频率响应的混叠失真，这是由于从s平面到z平面是多值的映射关系所造成的。
双线性变换法是使数字滤波器的频率响应与模拟滤波器的频率响应相似的一种变换方法。
为了克服多值映射这一缺点，我们首先把整个s平面压缩变换到某一中介的s1平面的一条横带里，再通过变换关系将此横带变换到整个z平面上去，这样就使得s平面与z平面是一一对应的关系，消除了多值变换性，也就消除了频谱混叠现象。
双线性法设计IIR数字滤波器的步骤：
1）将数字滤波器的频率指标{
}由Wk＝（2/T）*tan（
）,转换为模拟滤波器的频率指标{
}.
2）由模拟滤波器的指标设计H（s）.3）由H（s）转换为H（z）
4.22.FIR滤波器原理
FIR滤波器与IIR滤波器特点不同，IIR滤波器的相位是非线性的，若需线性相位则要采用全通网络进行相位校正。
而有限长单位冲激响应（FIR）数字滤波器就可以做成具有严格的线性相位，同时又可以具有任意的幅度特性。
由于FIR系统的冲激响应就是其系统函数各次项的系数，所以设计FIR滤波器的方法之一可以从时域出发，截取有限长的一段冲激响应作为H（z）的系数，冲激响应长度N就是系统函数H（z）的阶数。
只要N足够长，截取的方法合理，总能满足频域的要求。
这种时域设计、频域检验的方法一般要反复几个回合，不像IIRDF设计靠解析公式一次计算成功。
给出的理想滤波器频率响应是
，它是w的周期函数，周期
，因此可展开成傅氏级数
由傅立叶反变换导出，即
，再将
与窗函数
相乘就可以得到
，
。
、
的计算可采用傅氏变换的现成公式和程序，窗函数
也是现成的。
但整个设计过程不能一次完成，因为窗口类型和大小的选择没有解析公式可一次算，整个设计可用计算机编程来做。
窗函数的傅式变换W（ejω）的主瓣决定了H（ejω）过渡带宽。
W（ejω）的旁瓣大小和多少决定了H（ejω）在通带和阻带范围内波动幅度，常用的几种窗函数有：
矩形窗w（n）=RN（n）；
Hanning窗
；
Hamming窗
；
Blackmen窗
；
Kaiser窗
。
式中Io（x）为零阶贝塞尔函数。
五．设计步骤：
5.1录制女音：
利用MATLAB中的函数录制声音。
functionnvyin（）
fs=11025;%采样频率
str=['@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@'];
disp（str）;
disp（'开始录音'）;
str=['@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@'];
disp（str）;
y=wavrecord（3*fs,fs,'double'）;%录制声音3秒
str=['@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@'];
disp（str）;
disp（'录音结束'）;
str=['@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@'];
disp（str）;
str=['@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@'];
disp（str）;
disp（'播放录音'）;
str=['@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@'];
disp（str）;
wavplay（y,fs）;%播放录音
str=['@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@'];
disp（str）;
disp（'播放录音结束'）;
str=['@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@'];
disp（str）;
wavwrite（y,fs,'原女音'）;%声音的存储
5.2采样语音信号并画出时域波形和频谱图
读取语音信号，画出其时域波形和频谱图，与截取后的语音信号的时域波形和频谱图比较，观察其变化。
程序如下：
[x,fs,bits]=wavread（'女音.wav'）;%读取声音
N=length（x）;%计数读取信号的点数
t=（1:
N）/fs;%信号的时域采样点
f0=fs/N;%采样间隔
n=1:
N/2;%取信号的一半
figure
（1）;
subplot（2,2,1）;%把画图区域划分为2行2列，指定第一个图
plot（t,x）;%画出声音采样后的时域波形
title（'原女音信号的时域波形'）;%给图形加注标签说明
xlabel（'时间/t'）;
ylabel（'振幅/A'）;
grid;%添加网格
y=fft（x）;%对信号做N点FFT变换
k=（n-1）*f0;%频域采样点
subplot（2,2,3）;%把画图区域划分为2行2列，指定第三个图
plot（k,abs（y（n）））;%绘制原始语音信号的幅频响应图
title（'FFT变换后声音的频谱特性'）;%给图形加注标签说明
xlabel（'频率/Hz'）;
ylabel（'幅值/A'）;
grid;%添加网格
subplot（2,2,4）;%把画图区域划分为2行2列，指定第四个图
ify~=0%判断指数是否为0
plot（k,20*log10（abs（y（n））））;%画信号频谱的分贝图
end
xlabel（'Hz'）;
ylabel（'振幅/分贝'）;
title（'FFT变换后声音的频谱特性'）;
grid;%添加网格
%实际发出声音落后录制动作半拍的现象的解决
siz=wavread（'女音.wav','size'）;
x1=wavread（'女音.wav',[350032076]）;%截取语音信号
N=length（x1）;%计数读取信号的点数
t=（1:
N）/fs;%信号的时域采样点
f0=fs/N;%采样间隔
n=1:
N/2;%取信号的一半
figure
（2）;
subplot（2,2,1）;%把画图区域划分为2行2列，指定第一个图
plot（t,x1）;%画出声音采样后的时域波形
title（'截取后女音信号的时域波形'）;%给图形加注标签说明
xlabel（'时间/t'）;
ylabel（'振幅/A'）;
grid;%添加网格
y1=fft（x1）;%对信号做N点FFT变换
subplot（2,2,3）;%把画图区域划分为2行2列，指定第三个图
k=（n-1）*f0;%频域采样点
plot（k,abs（y（n）））;%绘制原始语音信号的幅频响应图
title（'FFT变换后声音的频谱特性'）;%给图形加注标签说明
xlabel（'频率/Hz'）;
ylabel（'幅值/A'）;
grid;%添加网格
subplot（2,2,4）;%把画图区域划分为1行2列，指定第二个图
ify1~=0%判断指数是否为0
plot（k,20*log10（abs（y1（n））））;%画信号频谱的分贝图
end
xlabel（'Hz'）;
ylabel（'振幅/分贝'）;
title（'FFT变换后声音的频谱特性'）;
grid;%添加网格
结果分析：
由原女音信号的时域波形可知录取开始时实际发出声音大概落后3500个采样点，我们把前3500点去除即可解决实际发出声音落后录制动作半拍的现象。
由原女音的的频谱图和截取后声音的频谱图可看出，对声音的截取并不会影响它们频谱分布。
5.3采用双线性变换法设计IIR滤波器：
人的声音频率一般在（1~~4）kHZ之间，则我们只需要设计一个带通滤波器即可滤去声音频带以外的无用噪声，得到比较清晰的声音。
根据声音的频谱图分析，设计一个带通滤波器性能指标如下：
fp1＝1000Hz，fp2＝3000Hz，fsc1＝500Hz，fsc2＝3500Hz，As＝100dB，Ap＝1dB，fs=10000
程序如下：
%iir带通的代码：
%w=2*pi*f/fs
Ap=1;%通带波纹系数
Az=100;%最小阻带衰减
wp=[0.20.6];%归一化通带数字截止频率
wz=[0.10.7];%归一化阻带数字截止频率
[N,wn]=cheb1ord（wp,wz,Ap,Az）;%估计契比雪夫I型滤波器阶数
[b,a]=cheby1（N,Ap,wn）;%N指定滤波器阶数，wn归一化截%止频率，Ap通带波动
[h,w]=freqz（b,a）;%求数字滤波器的复频率响应
figure
（1）;
subplot（2,1,1）;
plot（w/pi,abs（h））;%绘制数字滤波器的频谱图
grid;
xlabel（'\omega/\pi'）;
ylabel（'振幅（幅值）'）;
title（'契比雪夫Ⅰ型带通滤波器的幅频响应'）;
subplot（2,1,2）;
ifabs（h）~=0%判断指数是否为0
plot（w/pi,20*log10（abs（h）））;%绘制数字滤波器频谱的分贝图
end
grid;
xlabel（'\omega/\pi'）;
ylabel（'振幅（分贝）'）;
title（'契比雪夫Ⅰ型带通滤波器的幅频响应'）;
5.4窗函数法设计FFR滤波器
线性相位FIR滤波器通常采用窗函数法设计。
窗函数法设计FIR滤波器的基本思想是：
根据给定的滤波器技术指标，选择滤波器长度N和窗函数ω（n），使其具有最窄宽度的主瓣和最小的旁瓣。
其核心是从给定的频率特性，通过加窗确定有限长单位脉冲响应序列h（n）。
工程中常用的窗函数共有6种，即矩形窗、巴特利特（Bartlett）窗、汉宁（Hanning）窗、汉明（Hamming）窗、布莱克曼（Blackman）窗和凯泽（Kaiser）窗。
这次设计我采用的是布莱克曼来设计给定数字带通滤波器的参数如下：
wp1=0.3pi,wp2=0.6pi,wz1=0.2pi,wz2=0.7pi,Ap=1dB,Az=70dB
程序如下：
Ap=1;%通带波纹系数
Az=100;%最小阻带衰减
fs=10000;%采样频率
wp1=0.3*pi;
wp2=0.6*pi;
wz1=0.2*pi;
wz2=0.7*pi;
wc1=（wz1+wp1）/2;
wc2=（wz2+wp2）/2;
deltaW=min（（wp1-wz1）,（wz2-wp2））;%---取两个过渡带中的小者
N0=ceil（2*5.5*pi/deltaW）;%---查表7-3（P342）布拉克曼窗
N=N0+mod（N0+1,2）;%---确保N为奇数
hdWindow=ideallp（wc2,N）-ideallp（wc1,N）;%理想带通滤波器
wdWindow=blackman（N）;%布拉克曼窗
hr=wdWindow.*hdWindow';%点乘
n=0:
N-1;%阶数
subplot（2,2,1）;
stem（n,wdWindow）;%绘制布拉克曼窗时域波形
xlabel（'时间'）;
ylabel（'振幅'）;
title（'布拉克曼窗'）;
[H,W]=freqz（hr,1）;%求滤波器频率响应
subplot（2,2,3）;
plot（W/pi,abs（H））%绘制滤波器频域波形
xlabel（'\omega/\pi'）;
ylabel（'振幅'）;
title（'FIR带通滤波器幅频特性'）;
subplot（2,2,4）;
ifabs（H）~=0%判断指数是否为0
plot（W/pi,20*log10（abs（H）））;%画滤波器频谱的分贝图
end
xlabel（'\omega/'）;
ylabel（'振幅/分贝'）;
title（'FIR带通滤波器幅频特性'）;
grid;%添加网格
%---ideallp（）函数（非系统自有函数）在系统安装目录的WORK子目录ideallp.m
functionhd=ideallp（wc,N）;
%理想低通滤波器的脉冲响应子程序
%hd=点0到N-1之间的理想脉冲响应
%wc=截止频率（弧度）
%N=理想滤波器的长度
tao=（N-1）/2;%理想脉冲响应的对称中心位置
n=[0:
（N-1）];%设定脉冲响应长度
m=n-tao+eps;%加一个小数以避免零作除数
hd=sin（wc*m）./（pi*m）;%理想脉冲响应
5.5用IIR滤波器对信号进行滤波
用自己设计的IIR滤波器分别对采集的信号进行滤波，在Matlab中，IIR滤波器利用函数filter对信号进行滤波。
程序如下：
[x,fs,bits]=wavread（'女音.wav'）;
N=length（x）;%计数读取信号的点数
t=（1:
N）/fs;%信号的时域采样点
f0=fs/N;%采样间隔
n=1:
N/2;%取信号的一半
y=fft（x）;%对信号做N点FFT变换
k=（n-1）*f0;%频域采样点
subplot（2,1,1）;%把画图区域划分为2行1列，指定第一个图
plot（k,abs（y（n）））;%绘制原始语音信号的幅频响应图
title（'滤波前女音的频谱特性'）;%给图形加注标签说明
xlabel（'频率/Hz'）;
ylabel（'幅值/A'）;
grid;
%iir带通的代码：
Ap=1;%通带波纹系数
Az=100;%最小阻带衰减
wp=[0.20.6];%归一化通带数字截止频率
wz=[0.10.7];%归一化阻带数字截止频率
[N,wn]=cheb1ord（wp,wz,Ap,Az）;%估计契比雪夫I型滤波器阶数
[b,a]=cheby1（N,Ap,wn）;%N指定滤波器阶数，wn归一化截止频率，Ap通带波动
x1=filter（b,a,x）;%对声音滤波
wavplay（x1）
wavwrite（x1,'IIR滤波后女音.wav'）;
N=length（x1）;%计数读取信号的点数
t=（1:
N）/fs;%信号的时域采样点
f0=fs/N;%采样间隔
n=1:
N/2;%取信号的一半
y=fft（x1）;%对信号做N点FFT变换
k=（n-1）*f0;%频域采样点
subplot（2,1,2）;%把画图区域划分为2行1列，指定第一个图
plot（k,abs（y（n）））;%绘制原始语音信号的幅频响应图
title（'l滤波后女音的频谱特性'）;%给图形加注标签说明
xlabel（'频率/Hz'）;
ylabel（'幅值/A'）;
grid;
结果分析：
由上面滤波前后的频谱图可看出，滤波器滤除了小于1000Hz和大于3400Hz的频谱成分。
回放语音信号，由于低频和高频成分被滤除，声音变得较低沉。
5.6用FIR滤波器对信号进行滤波
用自己设计的FIR滤波器分别对采集的信号进行滤波，在Matlab中,FIR滤波器利用函数fftfilt对信号进行滤波
程序如下：
[x,fs,bits]=wavread（'女音.wav'）;
N=length（x）;%计数读取信号的点数
t=（1:
N）/fs;%信号的时域采样点
f0=fs/N;%采样间隔
n=1:
N/2;%取信号的一半
y=fft（x）;%对信号做N点FFT变换
k=（n-1）*f0;%频域采样点
subplot（2,1,1）;%把画图区域划分为2行1列，指定第一个图
plot（k,abs（y（n）））;%绘制原始语音信号的幅频响应图
title（'滤波前女音的频谱特性'）;%给图形加注标签说明
xlabel（'频率/\omega'）;
ylabel（'幅值/A'）;
grid;
%FIR带通滤波器代码
fs=10000;
wp1=0.3*pi;
wp2=0.6*pi;
wz1=0.2*pi;
wz2=0.7*pi;
wc1=（wz1+wp1）/2;
wc2=（wz2+wp2）/2;
deltaW=min（（wp1-wz1）,（wz2-wp2））;%---取两个过渡带中的小者
N0=ceil（2*5.5*pi/deltaW）;%---查表7-3（P342）布拉克曼窗
N=N0+mod（N0+1,2）;%---确保N为奇数
hdWindow=ideallp（wc2,N）-ideallp（wc1,N）;
wdWindow=blackman（N）;
hr=wdWindow.*hdWindow';
x1=fftfilt（hr,x）;%对声音滤波
wavplay（x1）
wavwrite（x1,'FIR滤波后女音.wav'）;
N=length（x1）;%计数读取信号的点数
t=（1:
N）/fs;%信号的时域采样点
f0=fs/N;%采样间隔
n=1:
N/2;%取信号的一半
y=fft（x1）;%对信号做N点FFT变换
k=（n-1）*f0;%频域采样点
subplot（2,1,2）;%把画图区域划分为2行1列，指定第一个图
plot（k,abs（y（n）））;%绘制原始语音信号的幅频响应图
title（'l滤波后女音的频谱特性'）;%给图形加注标签说明
xlabel（'频率/Hz'）;
ylabel（'幅值/A'）;
grid;
结果分析：
由上面滤波前后的频谱图可看出，滤波器滤除了小于1000Hz和大于3500Hz的频谱成分。
和用IIR滤波器滤波一样，回放语音信号，由于低频和高频成分被滤除，声音变得较低沉。
5.7男女声语音信号频谱特点分析
换一个男音录制同样一段语音内容，分析两段内容相同的语音信号频谱之间有什么特点。
程序如下：
[x,fs,bits]=wavread（'女音.wav'）;%读取声音
N=length（x）;%计数读取信号的点数
t=（1:
N）/fs;%信号的时域采样点
f0=fs/N;%采样间隔
n=1:
N/2;%取信号的一半
y=fft（x）;%对信号做N点FFT变换
k=（n-1）*f0;%频域采样点
subplot（2,1,1）;%把画图区域划分为2行1列，指定第一个图
plot（k,abs（y（n）））;%绘制原始语音信号的幅频响应图
title（'FFT变换后女音的频谱特性'）;%给图形加注标签说明
xlabel（'频率/\omega'）;
ylabel（'幅值/A'）;
grid;%添加网格
[x,fs,bits]=wavread（'明明.wav'）;%读取声音
N=length（x）;%计数读取信号的点数
t=（1:
N）/fs;%信号的时域采样点
f0=fs/N;%采样间隔
n=1:
N/2;%取信号的一半
y=fft（x）;%对信号做N点FFT变换
k=（n-1）*f0;%频域采样点
subplot（2,1,2）;%把画图区域划分为2行1列，指定第二个图
plot（k,abs（y（n）））;%绘制原始语音信号的幅频响应图
title（'FFT变换后男音的频谱特性'）;%给图形加注标签说明
xlabel（'频率/\omega'）;
ylabel（'幅值/A'）;
grid;%添加网格
axis（[060000300]）；%改变横纵坐标便于比较频谱图
结果分析：
通过比较上面女音频谱图和男音频谱图可知，男音的频谱集中在低频部分，高频成分底，谱线较平滑，声音听起来低沉。
5.8有背景噪声的信号分析
从硬盘中把一段噪声（频谱能量集中在某个小范围内）叠加到语音信号中，分析叠加前后信号频谱的变化，设计一个合适的滤波器，能够把该噪声滤除；
程序如下：
z=wavread（'女音.wav',[124000]）;%读取声音在1-24000之间
f=wavread（'noise.wav',[124000]）;
x=z+f;
wavplay（x）;
fs=11025;
N=length（x）;
f0=fs/N;%采样间隔
n=1:
N;%取信号的一半
y=fft（x,N）;%对信号做N点FFT变换
k=（n
升级会员 