数字信号处理课程设计说明书.docx
《数字信号处理课程设计说明书.docx》由会员分享，可在线阅读，更多相关《数字信号处理课程设计说明书.docx（28页珍藏版）》请在冰豆网上搜索。
数字信号处理课程设计说明书
广西科技大学
数字信号处理课程设计说明书
题目：
编程实现任意确定信号
的频谱分析算法（0.9）
系别：
计算机学院
专业：
通信工程
班级：
通信102班
学号：
201000402049
姓名：
黄绍耕
指导教师：
周坚和
日期：
2012年6月14日
目录
摘要………………………………………………………………….2
一、设计内容………………………………………………………2
二、设计原理………………………………………………………2
三、设计目的………………………………………………………4
四、实现过程………………………………………………………4
4.1.CEG和弦音音频文件的频谱分析…………………………4
4.2.对该信号频谱能量较集中的频带滤波……………………5
4.3.对滤波后的音频信号再滤出三个能量最集中的频簇……………9
4.4.重建信号与原信号的音频进行声音回放比较……………21
4.5.分析什么是和弦音…………………………………………22
五、课程设计总结…………………………………………………23
六、参考文献………………………………………………………23
摘要：
MATLAB主要是一种应用软件，可以运用高级语言编程思想，解决电子信息中的问题，可以做出仿真结果。
如信号处理，电路分析等。
本次课程设计主要是用MATLAB作为工具平台对给定的CEG和弦音音频文件进行分析处理，编程实现任意确定信号的频谱分析算法。
设计内容
（1）对给定的CEG和弦音音频文件取合适长度的采样记录点，然后进行频谱分析（信号的时域及幅频特性曲线要画出）。
（2）分析CEG和弦音频谱特点，对该信号频谱能量相对较为集中的频带（分低、中、高频）实现滤波（分别使用低通，带通及高通），显示滤波后信号的时域和频域曲线，并对滤波后的信号与原信号的音频进行声音回放比较。
（3）在低、中、高三个频带中，各滤出三个能量最集中的频簇，显示滤波后信号的时域和频域曲线。
（4）任意选择几个滤出的频带（或频簇）进行时域信号重建（合成），与原信号的音频进行声音回放比较。
讨论：
根据上述结果，分析什么是和弦音。
二、设计原理
1、采用双线性变换法设计滤波器，其原理如下：
S平面与z平面之间满足以下映射关系：
s平面的虚轴单值地映射于z平面的单位圆上，s平面的左半平面完全映射到z平面的单位圆内。
双线性变换不存在混叠问题。
双线性变换时一种非线性变换
，这种非线性引起的幅频特性畸变可通过预畸而得到校正。
IIR低通、高通、带通数字滤波器设计采用双线性原型变换公式：
变换类型
变换关系式
备 注
低通
高通
带通
为带通的上下
边带临界频率
可以利用上面提到的原理分别用双线性变化法设计以上3种滤波器，可以利用函数fir1设计FIR滤波器，可以利用函数butte，cheby1和ellip设计IIR滤波器；利用MATLAB中的函数freqz画出各滤波器的频率响应。
2、函数FFT用于序列快速傅立叶变换：
2.1函数的一种调用格式为：
y=fft（x）
其中，x是序列，y是序列的FFT，x可以为一向量或矩阵，若x为一向量，y是x的FFT。
且和x相同长度。
若x为一矩阵，则y是对矩阵的每一列向量进行FFT。
如果x长度是2的幂次方，函数fft执行高速基－2FFT算法；否则fft执行一种混合基的离散傅立叶变换算法，计算速度较慢。
2.2函数FFT的另一种调用格式为：
y=fft（x,N）式中，x，y意义同前，N为正整数。
函数执行N点的FFT。
若x为向量且长度小于N，则函数将x补零至长度N。
若向量x的长度大于N，则函数截短x使之长度为N。
若x为矩阵，按相同方法对x进行处理。
经函数fft求得的序列y一般是复序列，通常要求其幅值和相位。
MATLAB提供求复数的幅值和相位函数：
abs，angle，这些函数一般和FFT同时使用。
函数abs（x）用于计算复向量x的幅值，函数angle（x）用于计算复向量的相角，介于和之间，以弧度表示。
函数unwrap（p）用于展开弧度相位角p,当相位角绝对变化超过时，函数把它扩展至。
3、频率计算：
若N点序列x（n）（n=0,1,…,N-1）是在采样频率下获得的。
它的FFT也是N点序列，即X（k）（k=0,1,2,…,N-1），则第k点所对应实际频率值为f=k*fs/N.
三、设计目的：
综合运用本课程的理论知识进行频谱分析以及滤波器设计，通过理论推导得出相应结论，并利用MATLAB作为工具进行实现，从而复习巩固课堂所学的理论知识，提高对所学知识的综合应用能力，并从实践上初步实现对数字信号的处理。
四、实现过程
4.1.CEG和弦音音频文件的频谱分析
1、对给定的CEG和弦音音频文件取合适长度的采样记录点，然后进行频谱分析，并画出信号的时域及幅频特性曲线。
程序代码如下：
Fs=8000;%语音信号采样频率为8000
x=wavread（'C:
\MATLAB7\work\CEG和弦音.wav',[1008000]）;%把语音截取100到8000部分
sound（x）;%播放出语音
t=（0:
[8000-1000]）/Fs;%计算从1000到8000点的时间
y=fft（x,5000）;%对语音信号进行FFT运算
f=Fs*（0:
2499）/5000;
figure
（1）;
subplot（2,1,1）%在窗口中生成2行1列共2个子图，当前激活第1个图
plot（x）;%画出连续时域图形
title（'原和弦音信号的时域图'）;%标题为原始信号的时域图
xlabel（'时间'）;%定义图形横坐标为'时间’
ylabel（'幅值'）;%定义图形纵坐标为''幅值'
subplot（2,1,2）;%在窗口中生成2行1列共2个子图，当前激活第2个图
plot（f,abs（y（1:
2500）））;%画出幅值图
title（'原和弦音信号的时域图'）;%定义图形标题为'原和弦音信号的时域图'
xlabel（'频率'）;%定义图形横坐标为'频率’
ylabel（'幅值'）;%定义图形纵坐标为''幅值'
执行如下：
4.2.对该信号频谱能量较集中的频带滤波
2、对该信号频谱能量相对较为集中的频带（分低、中、高频）实现滤波（分别使用低通，带通及高通），同时显示滤波后信号的时域和频域曲线，并对滤波后的信号与原信号的音频进行声音回放比较。
程序代码如下：
（1）设计低通滤波器：
Fs=8000;
x=wavread（'C:
\MATLAB7\work\CEG和弦音.wav',[10008000]'）;%读出和弦音信号
t=（0:
[8000-1000]）/Fs;%计算从1000到8000点的时间
y=fft（x,5000）;%对和弦音信号进行FFT运算
f=Fs*（0:
2499）/5000;
Ws=2*1200*1/8000;%滤波器的阻带截止频率
Wp=2*1000*1/8000;%滤波器的通带截止频率
Rs=100;%定义通带最小阻带衰减
Rp=1;%定义通带波纹系数
[N,Wn]=cheb1ord（Wp,Ws,Rp,Rs）;%估计切贝雪夫I型滤波器阶数
[num,den]=cheby1（N,Rp,Wn,'low'）;%切贝雪夫I型低通滤波器系统函数；
[h,w]=freqz（num,den）;%数字低通滤波器的频率响应
subplot（3,1,1）;%在窗口中生成3行1列共3个子图，当前激活第1个子图
plot（w/pi,abs（h））;grid;%画出网格
xlabel（'\omega/\pi'）;%在xabel中显示归一化
/
ylabel（'振幅'）;%定义图形纵坐标为''振幅'
title（'契比雪夫Ⅰ型低通滤波器的幅频响应'）;%定义标题
f1=filter（num,den,x）;%滤波
y=fft（f1,8000）;%进行fft变换
subplot（3,1,2）;%在窗口中生成3行1列共3个子图，当前激活第2个子图
plot（f1）;%画出原语音信号经低通后图像
title（'低通滤波后的信号'）;%定义标题
xlabel（'时间'）;%定义图形横坐标为'时间’
ylabel（'幅值'）;%定义图形纵坐标为''幅值'
subplot（3,1,3）;%在窗口中生成3行1列共3个子图，当前激活第3个子图
plot（abs（y（1:
4000）））;%画出频谱图像
title（'低通后滤波信号频谱'）%定义标题
xlabel（'频率'）;%定义图形横坐标为'频率’
ylabel（'幅值'）;%定义图形纵坐标为''幅值'
wavwrite（f1,'低通.wav'）;%写出低通后信号
x1=wavread（'低通.wav'）;%读取低通后信号
sound（x1）;%播放低通后的语音
执行如下：
（2）设计高通滤波器：
程序代码如下：
Fs=8000;
x=wavread（'C:
\MATLAB7\work\CEG和弦音.wav',[10008000]）;%读出和弦音信号
t=（0:
[8000-1000]）/Fs;%计算从1000到8000点的时间
y=fft（x,5000）;%对和弦音信号进行FFT运算
f=Fs*（0:
2499）/5000;
Ws1=2*3000*1/8000;%滤波器的阻带截止频率
Wp1=2*3200*1/8000;%滤波器的通带截止频率
Rs1=100;%定义通带最小阻带衰减
Rp1=1;%定义通带波纹系数
[N1,Wn1]=cheb1ord（Wp1,Ws1,Rp1,Rs1）;%估计切贝雪夫Ⅰ型滤波器阶数
[num1,den1]=cheby1（N1,Rp1,Wn1,'high'）;%切贝雪夫I型高通滤波器系统函数；
[h1,w1]=freqz（num1,den1）;%计算幅频响应
subplot（3,1,1）;%在窗口中生成3行1列共3个子图，当前激活第1个子图
plot（w1/pi,abs（h1））;%画出幅值图
xlabel（'\omega/\pi'）;%在xabel中显示归一化
/
ylabel（'振幅'）;%定义图形纵坐标为''振幅'
title（'契比雪夫Ⅰ型高通滤波器的幅频响应'）;%定义标题
f2=filter（num1,den1,x）;%滤波
y1=fft（f2,8000）;%进行fft变换
subplot（3,1,2）;%在窗口中生成3行1列共3个子图，当前激活第2个子图
plot（f2）;%做原始语音信号的时域图形
title（'高通滤波后的信号'）;%定义标题
xlabel（'时间'）;%定义图形横坐标为'时间’
ylabel（'幅值'）;%定义图形纵坐标为''幅值'
subplot（3,1,3）;%在窗口中生成3行1列共3个子图，当前激活第3个子图
plot（abs（y1（1:
4000）））;%画出频谱图像
title（'高通滤波后的信号频谱'）;%定义标题
xlabel（'频率'）;%定义图形横坐标为'频率’
ylabel（'幅值'）;%定义图形纵坐标为'幅值'
wavwrite（f2,'高通.wav'）;%写出高通后信号
x2=wavread（'高通.wav'）;%读出高通后信号
sound（x2）;%播放高通后信号
执行如下：
（3）设计带通滤波器：
程序代码如下：
Fs=8000;
x=wavread（'C:
\MATLAB7\work\CEG和弦音.wav',[10008000]）;%读出和弦音信号
t=（0:
[8000-1000]）/Fs;%计算从1000到8000点的时间
y=fft（x,5000）;%对和弦音信号进行FFT运算
f=Fs*（0:
2499）/5000;
Wp2=[2*1200/Fs2*3000/Fs];%在1的dB衰减处的边带频率
Ws2=[2*1000/Fs2*3200/Fs];%在衰减为100dB处的边带频率
Rp2=1;%通带损耗不大于1dB
Rs2=100;%阻带衰减不小于100dB
[N2,Wn2]=cheb1ord（Wp2,Ws2,Rp2,Rs2）;%估计切贝雪夫滤波器阶数
[num2,den2]=cheby1（N2,Rp2,Wn2）;%切贝雪夫滤波器系统函数
[h2,w2]=freqz（num2,den2）;%计算频谱响应
figure（4）;%第四个图形
subplot（3,1,1）;%在窗口中生成3行1列共3个子图，当前激活第1个子图
plot（abs（h2））;grid;%画出幅值图
xlabel（'\omega/\pi'）;%在xabel中显示归一化
/
ylabel（'振幅'）;%定义图形纵坐标为''振幅'
title（'契比雪夫Ⅰ型带通滤波器的幅频响应'）;%定义标题
f3=filter（num2,den2,x）;%滤波
y3=fft（f3,8000）;%作fft变换
subplot（3,1,2）;%在窗口中生成3行1列共3个子图，当前激活第2个子图
plot（f3）;%做原始语音信号的时域图形
title（'带通滤波后的信号'）;%定义标题
xlabel（'时间'）;%定义图形横坐标为'时间’
ylabel（'幅值'）;%定义图形纵坐标为''幅值'
subplot（3,1,3）;%在窗口中生成3行1列共3个子图，当前激活第2个子图
plot（abs（y3（1:
4000）））;%画出频谱图像
title（'带通滤波后的信号频谱'）%定义标题
xlabel（'频率'）;%定义图形横坐标为'频率’
ylabel（'幅值'）;%定义图形纵坐标为''幅值'
wavwrite（f3,'带通.wav'）;%写出过滤后带通信号
x3=wavread（'带通.wav'）;%读出过滤后带通信号
sound（x3）;%播放带通后信号
执行如下：
4.3.对滤波后的音频信号再滤出三个能量最集中的频簇
3、在低、中、高三个频带中，各滤出三个能量最集中的频簇，显示滤波后信号的时域和频域曲线。
程序代码如下：
（1）低通能量最集中的频簇：
1、x1=wavread（'低通.wav'）;%读取低通后信号
sound（x1）;%播放低通后信号
Wp1=[2*500/Fs2*600/8000];
Ws1=[2*400/Fs2*800/8000];
Rp1=1;%通带损耗不大于1dB
Rs1=100;%阻带衰减不小于100dB
[N1,Wn1]=cheb1ord（Wp1,Ws1,Rp1,Rs1）;
[num1,den1]=cheby1（N1,Rp1,Wn1）;
[h1,w1]=freqz（num1,den1）;
figure（4）;%第四个图形
subplot（3,1,1）;%在窗口中生成3行1列共3个子图，当前激活第1个子图
plot（w1/pi,abs（h1））;grid;%打网格
xlabel（'\omega/\pi'）;%在xabel中显示归一化
/
ylabel（'振幅（幅值）'）;%定义图形纵坐标为''幅值'
title（'契比雪夫Ⅰ型低通滤波器的幅频响应'）;
f1=filter（num1,den1,x1）;%滤波
y1=fft（f1,8000）;%作fft变换
subplot（3,1,2）;%在窗口中生成3行1列共3个子图，当前激活第2个子图
plot（f1）;%做原始语音信号的时域图形
title（'低通能量集中后信号'）;
xlabel（'时间'）;%定义图形横坐标为'时间’
ylabel（'幅值'）;%定义图形纵坐标为''幅值'
subplot（3,1,3）;%在窗口中生成3行1列共3个子图，当前激活第3个子图
plot（abs（y1（1:
4000）））;%画出频谱图像
title（'低通的集中能量语音信号频谱'）%定义标题
xlabel（'频率'）;%定义图形横坐标为'频率’
ylabel（'幅值'）;%定义图形纵坐标为''幅值'
wavwrite（f1,'低通能量集中1.wav'）;%写出信号
x21=wavread（'低通能量集中1.wav'）;%读出信号
sound（x21）;%播放出信号
执行如下：
2、x2=wavread（'低通.wav'）;%读取低通后信号
sound（x2）;%播放低通后信号
Wp2=[2*700/80002*800/8000];
Ws2=[2*500/80002*1000/8000];
Rp2=1;
Rs2=100;
[N2,Wn2]=cheb1ord（Wp2,Ws2,Rp2,Rs2）;
[num2,den2]=cheby1（N2,Rp2,Wn2）;
[h2,w2]=freqz（num2,den2）;
figure（4）;
subplot（3,1,1）;%在窗口中生成3行1列共3个子图，当前激活第1个子图
plot（w2/pi,abs（h2））;grid;%打网格
xlabel（'\omega/\pi'）;%在xabel中显示归一化
/
ylabel（'振幅'）;%定义图形纵坐标为''振幅'
title（'契比雪夫Ⅰ型低通滤波器的幅频响应'）;
f2=filter（num2,den2,x2）;%滤波
y2=fft（f2,8000）;
subplot（3,1,2）;%在窗口中生成3行1列共3个子图，当前激活第2个子图
plot（f2）;
title（'低通能量集中后信号'）;
xlabel（'时间'）;%定义图形横坐标为'时间’
ylabel（'幅值'）;%定义图形纵坐标为''幅值'
subplot（3,1,3）;%在窗口中生成3行1列共3个子图，当前激活第3个子图
plot（abs（y2（1:
4000）））;%画出频谱图像
title（'低通的集中能量语音信号频谱'）
xlabel（'频率'）;%定义图形横坐标为'频率’
ylabel（'幅值'）;%定义图形纵坐标为''幅值'
wavwrite（f2,'低通能量集中2.wav'）;%写出信号
x31=wavread（'低通能量集中2.wav'）;%读出信号
sound（x31）;%播放出信号
执行如下：
3、x3=wavread（'低通.wav'）;%读取低通后信号
sound（x3）;%播放低通后信号
Wp3=[2*800/80002*900/8000];
Ws3=[2*600/80002*1200/8000];
Rp3=1;
Rs3=100;
[N3,Wn3]=cheb1ord（Wp3,Ws3,Rp3,Rs3）;
[num3,den3]=cheby1（N3,Rp3,Wn3）;
[h3,w3]=freqz（num3,den3）;
figure（4）;
subplot（3,1,1）;%在窗口中生成3行1列共3个子图，当前激活第1个子图
plot（w3/pi,abs（h3））;grid;%打网格
xlabel（'\omega/\pi'）;%在xabel中显示归一化
/
ylabel（'振幅'）;%定义图形纵坐标为''振幅
title（'契比雪夫Ⅰ型低通滤波器的幅频响应'）;
f3=filter（num3,den3,x3）;%滤波
y3=fft（f3,8000）;%做fft变换
subplot（3,1,2）;%在窗口中生成3行1列共3个子图，当前激活第2个子图
plot（f3）;
title（'低通能量集中后信号'）;
xlabel（'时间'）;%定义图形横坐标为'时间’
ylabel（'幅值'）;%定义图形纵坐标为''幅值
subplot（3,1,3）;%在窗口中生成3行1列共3个子图，当前激活第3个子图
plot（abs（y3（1:
4000）））;
title（'低通的集中能量语音信号频谱'）
xlabel（'频率'）;%定义图形横坐标为'频率’
ylabel（'幅值'）;%定义图形纵坐标为''幅值
wavwrite（f3,'低通能量集中3.wav'）;
x41=wavread（'低通能量集中3.wav'）;
sound（x41）;%播放出语音
执行如下：
（2）高通能量最集中的频簇：
1、x1=wavread（'高通.wav'）;%读出高通后信号
sound（x1）;%播放高通后信号
Wp1=[2*3300/80002*3500/8000];
Ws1=[2*3100/80002*3700/8000];
Rp1=1;
Rs1=100;
[N1,Wn1]=cheb1ord（Wp1,Ws1,Rp1,Rs1）;
[num1,den1]=cheby1（N1,Rp1,Wn1）;
[h1,w1]=freqz（num1,den1）;
figure（4）;
subplot（3,1,1）;%在窗口中生成3行1列共3个子图，当前激活第1个子图
plot（w1/pi,abs（h1））;grid;
xlabel（'\omega/\pi'）;%在xabel中显示归一化
/
ylabel（'振幅'）;%定义图形纵坐标为''幅值
title（'契比雪夫Ⅰ型高通滤波器的幅频响应'）;
f1=filter（num1,den1,x1）;%滤波
y1=fft（f1,8000）;%作fft变换
subplot（3,1,2）;%在窗口中生成3行1列共3个子图，当前激活第2个子图
plot（f1）;
title（'高通能量集中信号'）;%定义标题
xlabel（'时间'）;%定义图形横坐标为'时间’
ylabel（'幅值'）;%定义图形纵坐标为''幅值
subplot（3,1,3）;%在窗口中生成3行1列共3个子图，当前激活第3个子图
plot（abs（y1（1:
4000）））;
title（'高通能量集中的语音信号频谱'）%定义标题
xlabel（'频率'）;%定义图形横坐标为'频率’
ylabel（'幅值'）;%定义图形纵坐标为''幅值
wavwrite（f1,'高通后高通能量集中1.wav'）;
x22=wavread（'高通后高通能量集中1.wav'）;
sound（x22）;%播放出语音
执行如下：
2、x2=wavread（'高通.wav'）;%读出高通后信号
sound（x2）;%播放高通后信号
Wp2=[2*3000/80002*3200/8000];
Ws2=[2*2800/80002*3400/8000];
Rp2=1;
Rs2=100;
[N2,Wn2]=cheb1ord（Wp2,Ws2,Rp2,Rs2）;
[num2,den2]=cheby1（N2,Rp2,Wn2）;
[h2,w2]=freqz（num2,den2）;
figure（4）;
subplot（3,1,
升级会员 