编程实现任意确定信号的频谱分析算法.docx
《编程实现任意确定信号的频谱分析算法.docx》由会员分享，可在线阅读，更多相关《编程实现任意确定信号的频谱分析算法.docx（22页珍藏版）》请在冰豆网上搜索。
编程实现任意确定信号的频谱分析算法
IMBstandardizationoffice【IMB5AB-IMBK08-IMB2C】
编程实现任意确定信号的频谱分析算法
广西科技大学
数字信号处理课程设计说明书
题目：
编程实现任意确定信号的频谱分析算法
系别：
计算机工程学院
专业班级：
通信
学号：
学生姓名：
指导教师：
摘要…………………………………………………………………3
一、设计内容………………………………………………………3
二、设计原理………………………………………………………4
三、设计过程………………………………………………………7
和弦音音频文件的频谱分析………………………………7
2.对该信号频谱能量较集中的频带滤波…………………………9
3.比较滤波前后的音频文件………………………………………21
4.对滤波后的音频信号再滤出三个能量最集中的频簇…………21
5.重建信号与原信号的音频进行声音回放比较…………………33
6.分析什么是和弦音………………………………………………35
收获…………………………………………………………………36
参考文献……………………………………………………………36
摘要：
随着计算机和信息科学的飞速发展，信号处理逐渐发展成为一门独立的学科，成为信息科学的重要组成部分，在语音处理、雷达、图像处理、通信、生物医学工程等众多领域中得到广泛应应用。
Matlab语言是一种广泛应用于工程计算及数值分析领域的新型高级语言，Matlab功能强大、简单易学、编程效率高，深受广大科技工作者的喜爱。
特别是Matlab还具有信号分析工具箱，不需具备很强的编程能力，就可以很方便地进行语音信号分析、处理和设计。
数字信号处理课程在现代科学中具有很大重要性及自身特点，理解与掌握课程中的基本概念、基本原理、基本分析方法，对用Matlab进行数字信号处理课程设计的思路，具有很大帮助。
语音信号的处理与滤波的设计主要是用Matlab作为工具平台，设计中涉及到声音的录制、播放、存储和读取，语音信号的抽样、频谱分析，滤波器的设计及语音信号的滤波，通过数字信号处理课程的理论知识的综合运用。
从实践上初步实现对数字信号的处理。
关键词：
抽样频率；频谱分析；滤波器；窗函数
一、设计内容
（1）对给定的CEG和弦音音频文件取合适长度的采样记录点，然后进行频谱分析（信号的时域及幅频特性曲线要画出）。
（2）分析CEG和弦音频谱特点，对该信号频谱能量相对较为集中的频带（分低、中、高频）实现滤波（分别使用低通，带通及高通），显示滤波后信号的时域和频域曲线，并对滤波后的信号与原信号的音频进行声音回放比较。
（3）在低、中、高三个频带中，各滤出三个能量最集中的频簇，显示滤波后信号的时域和频域曲线。
（4）任意选择几个滤出的频带（或频簇）进行时域信号重建（合成），与原信号的音频进行声音回放比较。
讨论：
根据上述结果，分析什么是和弦音
二、设计原理
数字滤波器根据其冲激响应函数的时域特性，可分为两种，即无限长冲激响应（IIR）滤波器和有限长冲激响应（FIR）滤波器。
IIR滤波器的特征是，具有无限持续时间冲激响应。
种滤波器一般需要用递归模型来实现，因而有时也称之为递归滤波器。
FIR滤波器的冲激响应只能延续一定时间，在工程实际中可以采用递归的方式实现，也可以采用非递归的方式实现。
数字滤波器的设计方法有多种，如双线性变换法、窗函数设计法、插值逼近法和Chebyshev逼近法等等。
随着MATLAB软件尤其是MATLAB的信号处理工作箱的不断完善，不仅数字滤波器的计算机辅助设计有了可能，而且还可以使设计达到最优化.
线性相位FIR滤波器通常采用窗函数法设计。
窗函数法设计FIR滤波器的基本思想是：
根据给定的滤波器技术指标，选择滤波器长度N和窗函数ω（n），使其具有最窄宽度的主瓣和最小的旁瓣。
其核心是从给定的频率特性，通过加窗确定有限长单位脉冲响应序列h（n）。
工程中常用的窗函数共有6种，即矩形窗、巴特利特（Bartlett）窗、汉宁（Hanning）窗、汉明（Hamming）窗、布莱克曼（Blackman）窗和凯塞（Kaiser）。
FIR数字滤波器的设计原理：
采用窗口法
线性相位实系数FIR滤波器按其N值奇偶和h（n）的奇偶对称性分为四种：
1、h（n）为偶对称，N为奇数
H（ejω）的幅值关于ω=0，π，2π成偶对称。
2、h（n）为偶对称，N为偶数
H（ejω）的幅值关于ω=π成奇对称，不适合作高通。
3、h（n）为奇对称，N为奇数
H（ejω）的幅值关于ω=0，π，2π成奇对称，不适合作高通和低通。
4、h（n）为奇对称，N为偶数
H（ejω）ω=0、2π＝0，不适合作低通。
窗函数法设计线性相位FIR滤波器步骤:
确定数字滤波器的性能要求：
临界频率{ωk}，滤波器单位脉冲响应长度N；
根据性能要求，合理选择单位脉冲响应h（n）的奇偶对称性，从而确定理想频率响应Hd（ejω）的幅频特性和相频特性；
求理想单位脉冲响应hd（n），在实际计算中，可对Hd（ejω）按M（M远大于N）点等距离采样，并对其求IDFT得hM（n），用hM（n）代替hd（n）；
选择适当的窗函数w（n），根据h（n）=hd（n）w（n）求所需设计的FIR滤波器单位脉冲响应；
求H（ejω），分析其幅频特性，若不满足要求，可适当改变窗函数形式或长度N，重复上述设计过程，以得到满意的结果。
窗函数的傅式变换W（ejω）的主瓣决定了H（ejω）过渡带宽。
W（ejω）的旁瓣大小和多少决定了H（ejω）在通带和阻带范围内波动幅度，常用的几种窗函数有：
矩形窗w（n）=RN（n）；
Hanning窗
；
Hamming窗
；
Blackmen窗
；
Kaiser窗
。
式中Io（x）为零阶贝塞尔函数
三、设计过程
和弦音音频文件的频谱分析
对给定的CEG和弦音音频文件进行频谱分析，画出信号的时域及幅频特性曲线，程序如下：
[y,fs,bits]=wavread（''）;%采集已录声音信号
N=length（y）;%长度
t=（0:
length（y）-1）/fs;%计算时间
figure;
plot（t,y）;
xlabel（'时间'）;ylabel（'振幅'）;
title（'语音信号的时域图'）;
Y=fft（y）;%对信号进行FFT
df=fs/N;%采样间隔
n1=0:
N/2-1;
f=n1*df;%频带宽度
figure;
plot（f,abs（Y（1:
length（f））））;
xlabel（'频率'）;ylabel（'振幅'）;
title（'语音信号的频谱图'）;
sound（y）;%播放给定的CEG和谐音频文件
图形如下：
2.对该信号频谱能量较集中的频带滤波
对该信号频谱能量相对较为集中的频带（分低、中、高频）实现滤波（分别使用低通，带通及高通）
（1）FIR低通
Ws=2*1200*pi/8000;Wp=2*1000*pi/8000;As=100;Rp=1;%数字低通滤波器参数
Wc=（Wp+Ws）/2;%Wc
deltaW=Ws-Wp;%过度带ΔW取两个过渡带中的小者
N0=ceil*pi/deltaW）;N=N0+mod（N0+1,2）;%根据表格7-3求出响应的窗口的点数N，保证N为奇数
hd=ideallp（Wc,N）;%低通滤波器的hd（n）
%---ideallp（）函数（非系统自有函数）在系统安装目录的WORK子目录中，
wdWindow=kaiser（N,;%凯泽窗
hr=hd.*wdWindow';%---点乘
n=0:
N-1;
subplot（2,2,1）;
stem（n,wdWindow）;%画二维离散图形
xlabel（'时间'）;ylabel（'振幅'）;title（'凯泽窗'）;%x轴说明；y轴说明%标题
subplot（2,2,2）;
stem（n,hr）;%画二维离散图形
xlabel（'时间'）;ylabel（'振幅'）;%x轴说明；y轴说明
title（'凯泽窗的冲激响应'）;%标题
[H,W]=freqz（hr,1）;%求每个窗对应的频率响应
subplot（2,2,3）;plot（W/pi,abs（H））;%画二维连续图形
xlabel（'\omega/\pi'）;ylabel（'绝对振幅'）;%x轴说明；y轴说明
title（'滤波器低通幅频特性1'）;%标题
legend（'凯泽窗'）;
subplot（2,2,4）;
plot（W/pi,20*log10（abs（H）））;%画二维连续图形
xlabel（'\omega/\pi'）;ylabel（'振幅（分贝）'）;%x轴说明；y轴说明
title（'滤波器低通幅频特性2'）;%标题
legend（'凯泽窗'）;
[y,fs,bits]=wavread（''）;%读出给定的CEG和谐音频信号
N=length（y）;
t=（0:
length（y）-1）/fs;
df=fs/N;%采样间隔
n1=0:
N/2-1;
f=n1*df;%频带宽度
ly1=fftfilt（hr,y）;%信号送入滤波器滤波，ly1为输出
figure;%打开新的画图窗口
subplot（2,2,1）;
plot（t,y）;
xlabel（'时间'）;ylabel（'振幅'）;
title（'语音信号的时域图'）;
Y1=fft（y）;%滤波前的FFT变换
subplot（2,2,3）;
plot（f,abs（Y1（1:
length（f））））;%画滤波前频谱图
xlabel（'频率'）;ylabel（'振幅'）;
title（'语音信号的频谱图'）;
grid;
subplot（2,2,2）;
plot（t,ly1）;
xlabel（'时间'）;ylabel（'振幅'）;
title（'滤波后语音信号的时域图'）;
Ly1=fft（ly1）;%滤波后的FFT变换
subplot（2,2,4）;
plot（f,abs（Ly1（1:
length（f））））;%画滤波后频谱图
xlabel（'频率'）;ylabel（'振幅'）;
title（'滤波后语音信号的频谱图'）;
grid;
wavwrite（ly1,fs,'低通.wav'）;%把滤波后的信号写入文件
sound（ly1）;%回放低通滤波后的CEG和谐音频文件
（2）FIR带通
W1s=2*1000*pi/8000;W1p=2*1200*pi/8000;As=100;Rp=1;
W2s=2*3200*pi/8000;W2p=2*3000*pi/8000;%数字带通滤波器参数
W1c=（W1p+W1s）/2;
W2c=（W2s+W2p）/2;%W1c,W2c
deltaW=min（（W1p-W1s）,（W2s-W2p））;%过度带ΔW取两个过渡带中的小者
N0=ceil*pi/deltaW）;N=N0+mod（N0+1,2）;%根据表格7-3求出响应的窗口的点数N，保证N为奇数
hd=ideallp（W2c,N）-ideallp（W1c,N）;%带通滤波器的hd（n）
wdWindow=kaiser（N,;%凯泽窗
hr=hd.*wdWindow';%---点乘
n=0:
N-1;
subplot（2,2,1）;
stem（n,wdWindow）;%画二维离散图形
xlabel（'时间'）;ylabel（'振幅'）;%x轴说明；y轴说明
title（'凯泽窗'）;%标题
subplot（2,2,2）;
stem（n,hr）;%画二维离散图形
xlabel（'时间'）;ylabel（'振幅'）;%x轴说明；y轴说明
title（'凯泽窗的冲激响应'）;%标题
[H,W]=freqz（hr,1）;%求每个窗对应的频率响应
subplot（2,2,3）;plot（W/pi,abs（H））;%画二维连续图形
xlabel（'\omega/\pi'）;ylabel（'绝对振幅'）;%x轴，y轴说明
title（'滤波器带通幅频特性1'）;%标题
legend（'凯泽窗'）;
subplot（2,2,4）;
plot（W/pi,20*log10（abs（H）））;%画二维连续图形
xlabel（'\omega/\pi'）;ylabel（'振幅（分贝）'）;%x轴，y轴说明
title（'滤波器带通幅频特性2'）;%标题
legend（'凯泽窗'）;
[y,fs,bits]=wavread（''）;%读出给定的CEG音频文件
N=length（y）;
t=（0:
length（y）-1）/fs;
df=fs/N;%采样间隔
n1=0:
N/2-1;
f=n1*df;%频带宽度
ly2=fftfilt（hr,y）;%信号送入滤波器滤波，ly2为输出
figure;%打开新的画图窗口
subplot（2,2,1）;
plot（t,y）;
xlabel（'时间'）;ylabel（'振幅'）;
title（'语音信号的时域图'）;
Y2=fft（y）;%滤波前的FFT变换
subplot（2,2,3）;
plot（f,abs（Y2（1:
length（f））））;%画滤波前频谱图
xlabel（'频率'）;ylabel（'振幅'）;
title（'语音信号的频谱图'）;
grid;
subplot（2,2,2）;
plot（t,ly2）;
xlabel（'时间'）;
ylabel（'振幅'）;
title（'滤波后语音信号的时域图'）;
Ly2=fft（ly2）;%滤波后的FFT变换
subplot（2,2,4）;
plot（f,abs（Ly2（1:
length（f））））;%画滤波后频谱图
xlabel（'频率'）;ylabel（'振幅'）;
title（'滤波后语音信号的频谱图'）;
grid;
wavwrite（ly2,fs,'带通.wav'）;%把滤波后的信号写入文件
sound（ly2,fs）;%回放带通滤波后的CEG和谐音频文件
（3）FIR高通
Ws=2*3000*pi/8000;Wp=2*3200*pi/8000;As=100;Rp=1;%数字高通滤波器参数
Wc=（Wp+Ws）/2;%Wc
deltaW=Wp-Ws;%过度带ΔW取两个过渡带中的小者
N0=ceil*pi/deltaW）;N=N0+mod（N0+1,2）;%根据表格7-3求出响应的窗口的点数N，保证N为奇数
hd=ideallp（pi,N）-ideallp（Wc,N）;%高通滤波器的hd（n）
wdWindow=kaiser（N,;%凯泽窗
hr=hd.*wdWindow';%---点乘
n=0:
N-1;
subplot（2,2,1）;
stem（n,wdWindow）;%画二维离散图形
xlabel（'时间'）;ylabel（'振幅'）;%x轴,y轴说明
title（'凯泽窗'）;%标题
subplot（2,2,2）;
stem（n,hr）;%画二维离散图形
xlabel（'时间'）;ylabel（'振幅'）;%x轴.y轴说明
title（'凯泽窗的冲激响应'）;%标题
[H,W]=freqz（hr,1）;%求每个窗对应的频率响应
subplot（2,2,3）;plot（W/pi,abs（H））;%画二维连续图形
xlabel（'\omega/\pi'）;ylabel（'绝对振幅'）;%x轴,y轴说明
title（'滤波器高通幅频特性1'）;%标题
legend（'凯泽窗'）;
subplot（2,2,4）;plot（W/pi,20*log10（abs（H）））;%画二维连续图形
xlabel（'\omega/\pi'）;ylabel（'振幅（分贝）'）;%x轴,y轴说明
title（'滤波器高通幅频特性2'）;%标题
legend（'凯泽窗'）;
%在一个窗口同时画出滤波前后的波形及频谱
[y,fs,bits]=wavread（''）;%读出信号
N=length（y）;
t=（0:
length（y）-1）/fs;
df=fs/N;%采样间隔
n1=0:
N/2-1;
f=n1*df;%频带宽度
ly3=fftfilt（hr,y）;%信号送入滤波器滤波，ly3为输出
figure;%打开新的画图窗口
subplot（2,2,1）;
plot（t,y）;
xlabel（'时间'）;
ylabel（'振幅'）;
title（'语音信号的时域图'）;
Y3=fft（y）;%滤波前的FFT变换
subplot（2,2,3）;
plot（f,abs（Y3（1:
length（f））））;%画滤波前频谱图
xlabel（'频率'）;ylabel（'振幅'）;
title（'语音信号的频谱图'）;
grid;
subplot（2,2,2）;
plot（t,ly3）;
xlabel（'时间'）;ylabel（'振幅'）;
title（'滤波后语音信号的时域图'）;
Ly3=fft（ly3）;%滤波后的FFT变换
subplot（2,2,4）;
plot（f,abs（Ly3（1:
length（f））））;%画滤波后频谱图
xlabel（'频率'）;ylabel（'振幅'）;
title（'滤波后语音信号的频谱图'）;
grid;
wavwrite（ly3,fs,'高通.wav'）;%把滤波后的信号写入文件
sound（ly3,fs）;%回放高通滤波后的CEG音频文件
3.比较滤波前后的的音频文件
通过低、带、高通滤波后的音频文件与原音频文件有如下简单差异：
低通滤波后，CEG和谐音频文件的声音较原来变得低沉了。
带通滤波后，CEG和谐音频文件的声音没原来的那么尖锐。
高通滤波后，CEG和谐音频文件的声音变得尖锐清亮。
4对滤波后的音频信号再滤出三个能量最集中的频簇
实现低、带、高通滤波后，在三个频带中，再分别用带通滤波器滤出三个能量最集中的频簇，显示滤波后信号的时域和频域曲线。
（1）低通频簇
W1s=2*100*pi/8000;W1p=2*200*pi/8000;As=100;Rp=1;
W2s=2*900*pi/8000;W2p=2*800*pi/8000;%数字带通滤波器参数
W1c=（W1p+W1s）/2;
W2c=（W2s+W2p）/2;%W1c,W2c
deltaW=min（（W1p-W1s）,（W2s-W2p））;%过度带ΔW取两个过渡带中的小者
N0=ceil*pi/deltaW）;N=N0+mod（N0+1,2）;%根据表格7-3求出响应的窗口的点数N，保证N为奇数
hd=ideallp（W2c,N）-ideallp（W1c,N）;%带通滤波器的hd（n）
wdWindow=kaiser（N,;%凯泽窗
hr=hd.*wdWindow';%---点乘
n=0:
N-1;
subplot（2,2,1）;
stem（n,wdWindow）;%画二维离散图形
xlabel（'时间'）;ylabel（'振幅'）;%x轴，y轴说明
title（'凯泽窗'）;%标题
subplot（2,2,2）;
stem（n,hr）;%画二维离散图形
xlabel（'时间'）;ylabel（'振幅'）;%x轴，y轴说明
title（'凯泽窗的冲激响应'）;%标题
[H,W]=freqz（hr,1）;%求每个窗对应的频率响应
subplot（2,2,3）;
plot（W/pi,abs（H））;%画二维连续图形
xlabel（'\omega/\pi'）;ylabel（'绝对振幅'）;%x轴，y轴说明
title（'滤波器带通幅频特性1'）;%标题
legend（'凯泽窗'）;
subplot（2,2,4）;plot（W/pi,20*log10（abs（H）））;%画二维连续图形
xlabel（'\omega/\pi'）;ylabel（'振幅（分贝）'）;%x轴，y轴说明
title（'滤波器带通幅频特性2'）;%标题
legend（'凯泽窗'）;
[y1,fs,bits]=wavread（'低通.wav'）;%读出低通滤波后的信号
t=（0:
length（y1）-1）/fs;
df=fs/N;%采样间隔
n1=0:
N/2-1;
f=n1*df;%频带宽度
ly11=fftfilt（hr,y1）;%信号送入滤波器滤波，ly11为输出
figure;%打开新的画图窗口
subplot（2,2,1）;
plot（t,y1）;
xlabel（'时间'）;ylabel（'振幅'）;
title（'语音信号的时域图'）;
Y11=fft（y1）;%滤波前的FFT变换
subplot（2,2,3）;
plot（f,abs（Y11（1:
length（f））））;%画滤波前频谱图
xlabel（'频率'）;ylabel（'振幅'）;
title（'语音信号的频谱图'）;
grid;
subplot（2,2,2）;
plot（t,ly11）;
xlabel（'时间'）;ylabel（'振幅'）;
title（'滤波后语音信号的时域图'）;
Ly11=fft（ly11）;%滤波后的FFT变换
subplot（2,2,4）;
plot（f,abs（Ly11（1:
length（f））））;%画滤波后频谱图
xlabel（'频率'）;ylabel（'振幅'）;
title（'滤波后语音信号的频谱图'）;
grid;
wavwrite（ly11,fs,'低通'）;%把滤波后的信号写入文件
sound（ly11,fs）;%回放
（2）带通频簇
W1s=2*1100*pi/8000;W1p=2*1200*pi/8000;As=100;Rp=1;
W2s=2*2100*pi/8000;W2p=2*2000*pi/8000;%数字带通滤波器参数
W1c=（W1p+W1s）/2;
W2c=（W2s+W2p）/2;%W1c,W2c
deltaW=min（（W1p-W1s）,（W2s-W2p））;%过度带ΔW取两个过渡带中的小者
N0=ceil*pi/deltaW）;N=N0+mod（N0+1,2）;%根据表格7-3求出响应的窗口的点数N，保证N为奇数
hd=ideallp（W2c,N）-ideallp（W1c,N）;%带通滤波器的hd（n）
wdWindow=kaiser（N,;%凯泽窗
%---点乘
hr=hd.*wdWindow';
n=0:
N-1;
subplot（2,2,1）;
stem（n,wdWindow）;%画二维离散图形
xlabel（'时间'）;ylabel（'振幅'）;%x轴，y轴说明
title（'凯泽窗'）;%标题
subplot（2,2,2）;
stem（n,hr）;%画二维离散图形
xlabel（'时间'）;ylabel（'振幅'）;%x轴，y轴说明
title（'凯泽窗的冲激响应'）;%标题
[H,W]=freqz（hr,1）;%求每个窗对应的频率响应
subplot（2,2,3）;plot（W/pi,abs（H））;%画二维连续图形
xlabel（'\omega/\pi'）;ylabel（'绝对振幅'）;%x轴，y轴说明
title（'滤波器带通幅频特性1'）;%标题
l
升级会员 