完整word版综合训练.docx
《完整word版综合训练.docx》由会员分享，可在线阅读，更多相关《完整word版综合训练.docx（53页珍藏版）》请在冰点文库上搜索。
完整word版综合训练
数字信号处理课程研究报告
项目:
综合训练课题
班级:
测控142
姓名:
吉宇
学号:
160514205
学年:
2016—2017
常熟理工学院电气与自动化工程学院
课题一：
基于MATLAB的语音信号分析与处理
1、课题描述
录制一段自己的语音信号，并对录制的信号进行采样；画出采样后语音信号的时域波形和频谱波形；根据频率特征，分别设计IIR和FIR低通、高通、带通滤波器，并对声音信号进行滤波处理，回放声音信号，分析比较处理前后声音的变化。
二、课题分析
语音信号处理主要分成三个部分：
语音信号的录制与采样；
由于信号会有噪声，因此需要有相滤波器的设计；
最后对采样信号进行滤波。
三、课题设计
对原始声音信号采样并画出时域波形和频域波形程序设计如下：
原始声音信号采集
[y,fs,nbits]=wavread（'shengyin.wav'）;%读取声音
sound（y,fs,nbits）;%回放声音
N=length（y）;
n=0:
N-1;%长度
w=2*n*pi/N;
y1=fft（y）;%对原始信号做傅里叶变化
subplot（2,1,1）;
plot（n,y）;
title（'原始语音信号时域图'）;
xlabel（'时间t'）;
ylabel（'幅值y'）;
subplot（2,1,2）;
plot（w/pi,abs（y1））;
title（'原始语音信号频谱'）;
xlabel（'频率Hz'）;
ylabel（'幅度'）;
设计合适的滤波器前首先需要设计一个原型滤波器，下面分别设计了IIR和FIR低通，带通高通滤波器，设计程序过程如下：
原型滤波器设计
IIR低通滤波器：
Ft=8000;%模拟指标采样周期
Fp=1000;%通带截止频率
Fs=1200;%阻带截止频率
wp=2*pi*Fp/Ft;%数字指标
ws=2*pi*Fs/Ft;
ft=2*Ft*tan（wp/2）;%双线性变化
fs=2*Ft*tan（ws/2）;
[n11,wn11]=buttord（wp,ws,1,50,'s'）;
[b11,a11]=butter（n11,wn11,'s'）;
[num11,den11]=bilinear（b11,a11,0.5）;
[h,w]=freqz（num11,den11）;figure;
plot（w*8000*0.5/pi,abs（h））;
legend（'IIR低通滤波器'）;
grid;
IIR带通滤波器
Fp1=1200;
Fp2=3000;
Fs1=1000;
Fs2=3200;
Ft=8000;
wp1=tan（pi*Fp1/Ft）;
wp2=tan（pi*Fp2/Ft）;
ws1=tan（pi*Fs1/Ft）;
ws2=tan（pi*Fs2/Ft）;
w=wp1*wp2/ws2;
bw=wp2-wp1;
wp=1;
ws=（wp1*wp2-w.^2）/（bw*w）;
[n12,wn12]=buttord（wp,ws,1,50,'s'）;
[b12,a12]=butter（n12,wn12,'s'）;
[num2,den2]=lp2bp（b12,a12,sqrt（wp1*wp2）,bw）;
[num12,den12]=bilinear（num2,den2,0.5）;
[h,w]=freqz（num12,den12）;
figure;
plot（w*8000*0.5/pi,abs（h））;
axis（[0450001.5]）;
legend（'IIR带通滤波器','Location','NorthWest'）;
grid;
IIR高通滤波器：
Ft=8000;
Fp=4000;
Fs=3500;
wp1=tan（pi*Fp/Ft）;
ws1=tan（pi*Fs/Ft）;
wp=1;
ws=wp1*wp/ws1;
[n13,wn13]=cheb1ord（wp,ws,1,50,'s'）;
[b13,a13]=cheby1（n13,1,wn13,'s'）;
[num,den]=lp2hp（b13,a13,wn13）;
[num13,den13]=bilinear（num,den,0.5）;
[h,w]=freqz（num13,den13）;
figure;
plot（w*21000*0.5/pi,abs（h））;
legend（'IIR高通滤波器','Location','NorthWest'）;
axis（[01100001.5]）;
grid;
FIR低通滤波器：
Ft=8000;
Fp=1000;
Fs=1200;
wp=2*Fp/Ft;
ws=2*Fs/Ft;
rp=1;
rs=50;
p=1-10.^（-rp/20）;
s=10.^（-rs/20）;
fpts=[wpws];
mag=[10];
dev=[ps];
[n21,wn21,beta,ftype]=kaiserord（fpts,mag,dev）;
b21=fir1（n21,wn21,kaiser（n21+1,beta））;
[h,w]=freqz（b21,1）;
figure;
plot（w*8000*0.5/pi,abs（h））;
title（'FIR低通滤波器','fontweight','bold'）;
grid;
FIR带通滤波器
Fp1=1200;
Fp2=3000;
Fs1=1000;
Fs2=3200;
Ft=8000;
wp1=tan（pi*Fp1/Ft）;
wp2=tan（pi*Fp2/Ft）;
ws1=tan（pi*Fs1/Ft）;
ws2=tan（pi*Fs2/Ft）;
w=wp1*wp2/ws2;
bw=wp2-wp1;
wp=1;
ws=（wp*wp2-w.^2）/（bw*w）;
[n22,wn22]=buttord（wp,ws,1,50,'s'）;
[b22,a22]=butter（n22,wn22,'s'）;
[num2,den2]=lp2bp（b22,a22,sqrt（wp1*wp2）,bw）;
[num22,den22]=bilinear（num2,den2,0.5）;
[h,w]=freqz（num22,den22）;
figure;
plot（w*8000*0.5/pi,abs（h））;
axis（[0450001.5]）;
legend（'FIR带通滤波器','Locating','NorthWest'）;
grid;
FIR高通滤波器：
Ft=8001;
Fp=4000;
Fs=3500;
wp=2*Fp/Ft;
ws=2*Fs/Ft;
rp=1;
rs=50;
p=1-10.^（-rp/20）;
s=10.^（-rs/20）;
fpts=[ws,wp];
mag=[0,1];
dev=[p,s];
[n23,wn23,beta,ftype]=kaiserord（fpts,mag,dev）;
b23=fir1（n23,wn23,'high',kaiser（n23+1,beta））;
[h,w]=freqz（b23,1）;
figure;
plot（w*12000*0.5/pi,abs（h））;
title（'fir高通滤波器'）;
axis（[2500,5500,0,1.2]）;
grid;
对声音信号分别用所设计的滤波器进行滤波，分析滤波前后时域和频域波形的不同，设计程序如下：
对声音信号滤波处理后
IIR双线性低通滤波：
[y,fs,nbits]=wavread（'shengyin.wav'）;
n=length（y）;
S=fft（y）;
Ft=8000;
Fp=1000;
Fs=1200;
wp=2*pi*Fp/Ft;
ws=2*pi*Fs/Ft;
[n11,wn11]=buttord（wp,ws,1,50,'s'）;
[b11,a11]=butter（n11,wn11,'s'）;
[num11,den11]=bilinear（b11,a11,0.5）;
z11=filter（num11,den11,s）;
sound（z11,fs）;
m11=fft（z11）;
figure;
subplot（2,2,1）;
plot（abs（S）,'g'）;
title（'滤波前信号的频谱','fontweight','bold'）;
axis（[01500001000]）;
grid;
subplot（2,2,2）;
plot（abs（m11）,'r'）;
title（'滤波后信号的频谱','fontweight','bold'）;
axis（[01500001000]）;
grid;
subplot（2,2,3）;
plot（s）;
title（'滤波前信号的波形','fontweight','bold'）;
axis（[6700087000-0.50.5]）;
grid;
subplot（2,2,4）;
plot（z11）;
title（'滤波后信号的波形','fontweight','bold'）;
axis（[6700087000-0.50.5]）;
grid;
IIR双线性带通滤波：
[y,fs,nbits]=wavread（'shengyin.wav'）;
n=length（y）;
S=fft（y）;
Ft=20000;
Fp=700;
Fs=1400;
wp=2*Fp/Ft;
ws=2*Fs/Ft;
rp=1;
rs=50;
p=1-10.^（-rp/20）;
q=10.^（-rs/20）;
fpts=[wpws];
mag=[10];
dev=[pq];
[n21,wn21,beta,ftype]=kaiserord（fpts,mag,dev）;
b21=fir1（n21,wn21,kaiser（n21+1,beta））;
z21=fftfilt（b21,s）;
sound（z21,fs）;
m21=fft（z21）;
figure（4）;
subplot（2,2,1）;
plot（abs（S）,'g'）;
title（'滤波前信号的频谱','fontweight','bold'）;
axis（[01500001000]）;
grid;
subplot（2,2,2）;
plot（abs（m11）,'r'）;
title（'滤波后信号的频谱','fontweight','bold'）;
axis（[01500001000]）;
grid;
subplot（2,2,3）;
plot（s）;
title（'滤波前信号的波形','fontweight','bold'）;
axis（[6700087000-0.50.5]）;
grid;
subplot（2,2,4）;
plot（z11）;
title（'滤波后信号的波形','fontweight','bold'）;
axis（[6700087000-0.50.5]）;
grid;
IIR双线性高通滤波：
[y,fs,nbits]=wavread（'shengyin.wav'）;
n=length（y）;
S=fft（y）;
Fp1=1400;
Fs1=700;
Ft=10000;
wp1=tan（pi*Fp1/Ft）;
ws1=tan（pi*Fs1/Ft）;
wp=1;
ws=wp*wp/ws1;
[n13,wn13]=cheb1ord（wp,ws,1,50,'s'）;
[b13,a13]=cheby1（n13,1,wn13,'s'）;
[num,den]=lp2hp（b13,a13,wn13）;
[num13,den13]=bilinear（num,den,0.5）;
z13=filter（num13,den13,s）;
sound（z13,fs）;
m13=fft（z13）;
figure;
subplot（2,2,1）;
plot（abs（S）,'g'）;
title（'滤波前信号的频谱','fontweight','bold'）;
axis（[01500001000]）;
grid;
subplot（2,2,2）;
plot（abs（m11）,'r'）;
title（'滤波后信号的频谱','fontweight','bold'）;
axis（[01500001000]）;
grid;
subplot（2,2,3）;
plot（s）;
title（'滤波前信号的波形','fontweight','bold'）;
axis（[6700087000-0.50.5]）;
grid;
subplot（2,2,4）;
plot（z11）;
title（'滤波后信号的波形','fontweight','bold'）;
axis（[6700087000-0.50.5]）;
grid;
FIR窗函数低通滤波：
[y,fs,nbits]=wavread（'shengyin.wav'）;
n=length（y）
S=fft（y）;;
Ft=10000;
Fp=700;
Fs=1400;
wp=2*Fp/Ft;
ws=2*Fs/Ft;
rp=1;
rs=50;
p=1-10.^（-rp/20）;
q=10.^（-rs/20）;
fpts=[wpws];
mag=[10];
dev=[pq];
[n21,wn21,beta,ftype]=kaiserord（fpts,mag,dev）;
b21=fir1（n21,wn21,kaiser（n21+1,beta））;
z21=fftfilt（b21,s）;
sound（z21,fs）;
m21=fft（z21）;
figure（4）;
subplot（2,2,1）;
plot（abs（S）,'g'）;
title（'滤波前的信号频谱','fontweight','bold'）;
axis（[01500001000]）;
grid;
subplot（2,2,2）;
plot（abs（m21）,'r'）;
title（'滤波后的信号频谱','fontweight','bold'）;
axis（[01500001000]）;
grid;
subplot（2,2,3）;
plot（s）;
title（'滤波前的信号波形','fontweight','bold'）;
axis（[6700087000-0.50.5]）;
grid;
subplot（2,2,4）;
plot（z21）;
title（'滤波后的信号波形','fontweight','bold'）;
axis（[6700087000-0.50.5]）;
grid;
FIR窗函数带通滤波：
[y,fs,nbits]=wavread（'shengyin.wav'）;
n=length（y）
S=fft（y）;
Fp1=1200;
Fp2=3000;
Fs1=1000;
Fs2=3200;
Ft=2200;
wp1=tan（pi*Fp1/Ft）;
wp2=tan（pi*Fp2/Ft）;
ws11=tan（pi*Fs1/Ft）;
ws2=tan（pi*Fs2/Ft）;
w=wp1*wp2/ws2;
bw=wp2-wp1;
wp=1;
ws=（wp*wp2-w.^2）/（bw*w）;
[n22,wn22]=buttord（wp,ws,1,50,'s'）;
[b22,a22]=butter（n22,wn22,'s'）;
z22=fftfilt（b22,s）;
sound（z22,fs）;
m22=fft（z22）;
figure;
subplot（2,2,1）;
plot（abs（S）,'g'）;
title（'滤波前的信号频谱','fontweight','bold'）;
axis（[01500001000]）;
grid;
subplot（2,2,2）;
plot（abs（m22）,'r'）;
title（'滤波后的信号频谱','fontweight','bold'）;
axis（[01500001000]）;
grid;
subplot（2,2,3）;plot（s）;
title（'滤波前的信号波形','fontweight','bold'）;
axis（[67008700-0.50.5]）;
grid;
subplot（2,2,4）;
plot（z22）;
title（'滤波后的信号波形','fontweight','bold'）;
axis（[6700087000-0.50.5]）;
grid;
FIR窗函数高通滤波：
[y,fs,nbits]=wavread（'shengyin.wav'）;
n=length（y）
S=fft（y）;;
Ft=10000;
Fp=1400;
Fs=700;
wp=2*Fp/Ft;
ws=2*Fs/Ft;
rp=1;
rs=50;
p=1-10.^（-rp/20）;
q=10.^（-rs/20）;
fpts=[wswp];
mag=[10];
dev=[pq];
[n23,wn23,beta,ftype]=kaiserord（fpts,mag,dev）;
b23=fir1（n23,wn23,'high',kaiser（n23+1,beta））;
z23=fftfilt（b23,s）;
sound（z23,fs）;
m23=fft（z23）;
figure;
subplot（2,2,1）;
plot（abs（S）,'g'）;
title（'滤波前的信号频谱','fontweight','bold'）;
axis（[01500001000]）;
grid;
subplot（2,2,2）;
plot（abs（m23）,'r'）;
title（'滤波后的信号频谱','fontweight','bold'）;
axis（[01500001000]）;
grid;
subplot（2,2,3）;
plot（s）;
title（'滤波前的信号波形','fontweight','bold'）;
axis（[6700087000-0.50.5]）;
grid;
subplot（2,2,4）;
plot（z23）;
title（'滤波后的信号波形','fontweight','bold'）;
axis（[6700087000-0.50.5]）;
grid;
课题总结：
本次的数字信号处理综合实验的题目是应用Matlab对语音信号进行频谱分析及滤波，首先通过网络和书籍查找有关本次综合实验的材料，编写相关程序，并通过Matlab软件运行得到相关波形频谱图。
实验中利用双线性变换法设计IIR数字滤波器，利用窗函数设计FIR数字滤波器，可以是低通、高通和带通滤波器的设计，通过设计的滤波器对语音信号进行滤波，再对得出的频谱图进行分析。
在实验中得到一些困难，在设计数字滤波器的时候，通带和阻带频率的选取需要满足低通的要求，以及带通允许的最大衰减和阻带应达到最小的衰减。
课题二：
基于MATLAB的谐波分析FFT
本文应用MATLAB来验证定理：
方波可用相应频率的基波及其奇次谐波合成。
一、 公式分析及计算
1.1傅里叶变换的原理
任何具有性质周期为T的波函数（）ft都可以表示为三角函数所构成的级数之和，即：
（1） 其中：
t为时间，ω为角频率。
Ω=2Π/T（T为周期），第一项1/2（a0）为直流分量。
图1 方波
所谓周期性函数的傅里叶变换（Fourier transform）就是将周期性函数张凯成直流分量，基波和所有n次谐波的叠加。
图1所示的方波可以写成函数形式：
f（t）=h,（0-h,（-T/2≤t<0）
在这里，h为常数2。
很明显，此方波为奇函数，并且它没有常数项，同时，它是一个周期为T的函数，所以我们可以用傅里叶级数来表示这个函数。
我们把它展开，可以得到：
1.2傅里叶变换的证明
下面，我们要从数学角度来证明为什么公式（3）能成立。
由于这是一个奇函数,常数项0a可以用积分函数计算出来：
所以其常数项不存在，即a0=0，下面开始计算an与bn:
1.3 周期信号的分解
周期信号是定义在（-∞,+ ∞）区间，每隔一定的时间T，按相同规律重复变化的信号，它可表示为:
f（t）=f（t+mT）
式中m为任意整数。
时间T称为该信号的重复周期，简称为周期。
需要指出的是，只有当周期信号满足狄里赫利条件时，才能展开为傅里叶级数。
狄里赫利条件是：
1）函数在任意有限区间内连续，或只有有限个第一类间断点（当t从左或右趋于这个间断点时，函数有有限的左极限和右极限） 2）在一周期内，函数有有限个极大值或极小值。
上式中系数an,bn称为傅里叶系数。
为简便式，
积分区间（t0, t0+T）取为（-2T，2T）或（0，T）。
考虑到正、余弦函数的正交条件，可得傅里叶系数。
以下是公式（2-2），公式（2-3）
式中T为函数f（t）的周期，Ω=2Π/T为角频率，由上述两式，傅里叶系数an和bn都是n的函数，其中an是n的偶函数，即a（-n）= an；而bn是n的奇函数，既有b（-n）=- bn
将式（2-1）中同频率项合并，可写成如下形式:
如将式（2-4）的形式化为（2-1）的形式，他们系数之间的关系为:
式（2-4）表明，任何满足狄里赫利条件的周期函数可分解为直流和许多余弦（或正弦）分量。
其中第一项A/2是常数项，它是周期信号中所包涵的直流分量；式中第二项A1cos（1Ω+φ1）称为基波或一次谐波，它的角频率与原周期信号相同，A1是基波振幅，1是基波初相角；式中第三项A2cos（2Ω+φ1）称为二次谐波，它的频率是基波频率的两倍，A2是二次谐波振幅，2是其初相角。
以此类推，还有三次、四次、„„谐波。
一般而言，Ancos（nΩ+φ1）称为n次谐波，An是n次谐波的振幅，n是其初相角。
式（2-4）表明，周期函数可以分解为各谐波分量.
1.4 方波的分解
设方波信号f（t）的周期为T，宽度为2T，将其展开为傅里叶级数 由式（2-2）和（2-3）可得
将它们代入到式（2-1），得到信号的傅里叶级数展开式为
它只含一、三、五„奇次谐波分量。
下图中画出了一个周期的方波组成情况，由图可见，当它包含的谐波分量愈多时，波形就愈接近原来的方波信号f（t）（图中虚线所示），其均方误差愈小，还可以看出，频率较低的谐波，其振幅较大，他们组成方波的主体，而频率较高的高次谐波振幅较小，它们主要影响波形的细节，波形中所包含的高次谐
升级会员 