基于MATLAB的频谱分析与信号去噪讲解.docx
《基于MATLAB的频谱分析与信号去噪讲解.docx》由会员分享，可在线阅读，更多相关《基于MATLAB的频谱分析与信号去噪讲解.docx（29页珍藏版）》请在冰豆网上搜索。
基于MATLAB的频谱分析与信号去噪讲解
学生毕业设计报告
基于MATLAB的频谱分析与信号去噪
作 者
系（院）
专 业
年 级
学 号
指导教师
日 期
学生诚信承诺书
本人郑重承诺：
所呈交的设计报告是我个人在导师指导下进行的研究工作及取得的研究成果。
尽我所知，除了文中特别加以标注和致谢的地方外，报告中不包含其他人已经发表或撰写的研究成果，也不包含为获得安阳师范学院或其他教育机构的学位或证书所使用过的材料。
与我一同工作的同志对本研究所做的任何贡献均已在报告中作了明确的说明并表示了谢意。
签名：
日期：
报告使用授权说明
本人完全了解有关保留、使用学位报告的规定，即：
学校有权保留送交报告的复印件，允许报告被查阅和借阅；学校可以公布报告的全部或部分内容，可以采用影印、缩印或其他复制手段保存报告。
签名：
导师签名：
日期：
基于MATLAB的频谱分析与信号去噪
摘要：
本课题是基于数字信号处理的理论知识对语音信号、图像信号等的频谱分析以及对加噪声语音信号进行时域、频域分析和滤波设计。
然后利用MATLAB软件进行编程，调试并完善程序，最终在计算机上得以实现。
滤波器设计在数字信号处理中占有极其重要的地位，数字滤波器又有FIR数字滤波器和IIR滤波器两种。
利用MATLAB信号处理工具箱可以快速有效地设计各种数字滤波器。
在设计实现的过程中，使用窗函数法来设计FIR数字滤波器，用巴特沃斯、切比雪夫和双线性变法设计IIR数字滤波器，并利用MATLAB作为辅助工具完成设计中的计算与图形的绘制。
通过对对所设计滤波器的仿真和频率特性分析，可知利用MATLAB信号处理工具箱可以有效快捷地设计FIR和IIR数字滤波器，过程简单方便，结果的各项性能指标均符合指定要求。
关键词频域分析，滤波器，MATLAB
1引言
1.1课题背景
随着信息时代和数字世界的到来，数字信号处理己成为当今一门极其重要的学科和技术领域，数字信号处理在通信、语音、图像、自动控制、医疗和家用电器等众多领域得到了广泛的应用。
任意一个信号都具有时域与频域特性，信号的频谱完全代表了信号，因而研究信号的频谱就等于研究信号本身。
通常从频域角度对信号进行分析与处理，容易对信号的特性获得深入的了解。
因此，信号的频谱分析是数字信号处理技术中的一种较为重要的工具。
数字滤波器,是数字信号处理中及其重要的一部分。
随着数字技术的发展，受到人们越来越多的重视。
数字滤波器可以通过数值运算实现滤波，所以数字滤波器处理精度高、稳定、体积小、重量轻、灵活、不存在阻抗匹配问题，可以实现模拟滤波器无法实现的特殊功能。
数字滤波器种类很多，根据其实现的网络结构或者其冲激响应函数的时域特性，可分为两种，即有限冲激响应（FIR，FiniteImpulseResponse）滤波器和无限冲激响应（IIR，InfiniteImpulseResponse）滤波器。
在工程领域中，MATLAB是一种倍受程序开发人员青睐的语言，对于一些需要做大量数据运算处理的复杂应用以及某些复杂的频谱分析算法MATLAB显得游刃有余。
1.2研究意义
信号处理几乎涉及到所有的工程技术领域，而频谱分析又是信号处理中一个非常重要的分析手段。
一般的频谱分析都依靠传统频谱分析仪来完成，价格昂贵，体积庞大，不便于工程技术人员的携带。
虚拟频谱分析仪改变了原有频谱分析仪的整体设计思路，用软件代替了硬件，使工程技术人员可以用一部笔记本电脑到现场就可轻松完成信号的采集、处理及频谱分析。
信号去噪被用于从一堆波音资料中提取有用信息，去除干扰，提高波音资料信噪比。
为了提高信噪比，人们根据信号和噪声的各种特征差异，设计了许多去噪方法，并在应用中并取得了很好的成果。
信号去噪的很多方法都是利用短时傅立叶变换来滤波去噪，但是短时傅立叶变换不能同时兼顾时间分辨率和频率分辨率。
用不同种滤波器滤波也是一种有效可行的方法。
1.3本文研究内容
信号的频谱分析就是利用傅里叶分析的方法，求出与时域描述相对应的频域描述,从中找出信号频谱的变化规律，以达到特征提取的目的。
不同信号的傅里叶分析理论与方法，在有关专业书中都有介绍。
但实际的待分析信号一般没有解析式，直接利用公式进行傅里叶分析非常困难.。
DFT是一种时域和频域均离散化的傅里叶变换，适合数值计算且有快速算法，是分析信号的有力工具。
DFT及FFT是数字信号处理的重要内容。
DFT是FFT的基础，FFT是DFT的快速算法，在MATLAB中可以利用函数FFT来计算序列的离散傅里叶变换DFT。
基于此首先介绍了MATLAB软件；其次给出了基于MATLAB软件的DFT和FFT频谱分析的方法，利用MATLAB软件方法，使得设计方便、快捷，大大减轻了工作量；再而我们提取一段语音信号，加噪，再通过设计的不同滤波器进行滤波后的频谱分析。
本文将重点介绍基于MATLAB的频谱分析设计，包括：
（1）音频信号频谱分析；
（2）图像信号频谱分析；
（3）离散信号/序列频谱分析；
（4）语音信号提取，分析和加噪；
（5）设计IIR和FIR的各种滤波器；
（6）用设计的滤波器进行滤波；
（7）分析滤波前后信号特征。
2频谱分析技术及MATLAB简介
2.1时域抽样定理
时域抽样定理给出了连续信号抽样过程中信号不失真的约束条件：
对于基带信号，信号抽样频率
大于等于2倍的信号最高频率
，即
。
时域抽样是把连续信号
变成适于数字系统处理的离散信号X[k]。
对连续信号
以间隔T抽样，则可得到的离散序列为
，如图2-1所示。
图2-1连续信号抽样的离散序列
若
，则信号
与X[k]的频谱之间存在：
其中，
的频谱为
，X[k]的频谱为
。
可见，信号时域抽样导致信号频谱的周期化。
（rad/s）为抽样角频率，
为抽样频率。
数字角频率Ω与模拟角频率ω的关系为：
Ω=ωT。
2.2离散傅立叶变换（DFT）
有限长序列
的离散傅立叶变换（DFT）为
逆变换为
2.3快速傅立叶变换（FFT）
在各种信号序列中，有限长序列占重要地位。
对有限长序列可以利用离散傅立叶变换（DFT）进行分析。
DFT不但可以很好的反映序列的频谱特性，而且易于用快速算法（FFT）在计算机上进行分析。
有限长序列的DFT是其z变换在单位圆上的等距离采样，或者说是序列傅立叶的等距离采样，因此可以用于序列的谱分析。
FFT是DFT的一种快速算法，它是对变换式进行一次次分解，使其成为若干小数据点的组合，从而减少运算量。
MATLAB为计算数据的离散快速傅立叶变换，提供了一系列丰富的数学函数，主要有fft、ifft、fft2、ifft2,fftn、ifftn和fftshift、ifftshift等。
当所处理的数据的长度为2的幂次时，采用基-2算法进行计算，计算速度会显著增加。
所以，要尽可能使所要处理的数据长度为2的幂次或者用添零的方式来添补数据使之成为2的幂次。
fft函数调用方式：
（1）Y=fft（X）
（2）Y=fft（X,N）
（3）Y＝fft（X,[],dim）或Y＝fft（X,N,dim）。
函数ifft的参数应用与函数fft完全相同。
2.4频谱分析原理
时域分析只能反映信号的幅值随时间的变化情况，除单频率分量的简单波形外，很难明确提示信号的频率组成和各频率分量大小，而频谱分析能很好的解决此问题。
由于从频域能获得的主要是频率信息，所以本节主要介绍频率（周期）的估计与频谱图的生成。
2.4.1频率和周期的估计
对于Y（kΔf），如果当kΔf = f时，Y（kΔf）取最大值，则f为频率的估计值，由于采样间隔的误差，f也存在误差，其误差最大为Δf /2。
周期T=1/f。
从原理上可以看出，如果在标准信号中混有噪声，用上述方法仍能够精确地估计出原标准信号的频率和周期。
2.4.2频谱图
为了直观地表示信号的频率特性，工程上常常将Fourier变换的结果用图形的方式表示，即频谱图。
以频率f为横坐标，|Y（f）|为纵坐标，可以得到幅值谱；
以频率f为横坐标，Arg Y（f）为纵坐标，可以得到相位谱；
以频率f为横坐标，Re Y（f）为纵坐标，可以得到实频谱；
以频率f为横坐标，Im Y（f）为纵坐标，可以得到虚频谱。
根据采样定理，只有频率不超过Fs/2的信号才能被正确采集，即Fourier变换的结果中频率大于Fs/2的部分是不正确的部分，故不在频谱图中显示。
即横坐标f ∈[0, Fs/2]
2.5MATLAB简介
2.5.1MATLAB软件的发展
MATLAB软件是由美国Mathworks公司推出的用于数值计算和图形处理的科学计算系统环境。
MATLAB是英文MATrixLABoratory（矩阵实验室）的缩写。
它的第一版（DOS版本1.0）发行于1984年，经过10余年的不断改进，现今已推出它的Windows98/NT版本（6.1版）。
新的版本集中了日常数字处理中的各种功能，包括高效的数值计算、矩阵运算、信号处理和图形生成等功能。
在MATLAB环境下，用户可以集成地进行程序设计、数值计算、图形绘制、输入输出、文件管理等各项操作。
MATLAB提供了一个人机交互的数学系统环境，该系统的基本数据结构是矩阵，在生成矩阵对象时，不要求作明确的维数说明。
与利用C语言或FORTRAN语言作数值计算的程序设计相比，利用MATLAB可以节省大量的编程时间。
在美国的一些大学里，MATLAB正在成为对数值计算、算法预设计与验证，以及一些特殊的短阵计算应用，如自动控制理论、统计、数字信号处理（时间序列分拆）等。
MATLAB系统最初是由CieveMoler用FORTRAN语言设计的，有关矩阵的算法来自LINPACK和EISPACK课题的研究成果；现在的MATLAB程序是MathWorks公司用C语言开发的，第一版由steveBangert主持开发编译解释程序，SteveKleiman完成图形功能的设计，JohnLittle和CleveMoler主持开发各类数学分析的子模块，撰写用户指南和大部分M文件。
自从第1版发行以来，已有众多的科技工作者加入到MATLAB的开发队伍中，并为形成今天的MATLAB系统做出了巨大的贡献，MATLAB以商品形式出现后，仅短短几年，就以其良好的开放性和运行可靠性，使原先控制领域里的封闭式软件包（如英国的UMIST，瑞典的LUND，德国的KEDDC）纷纷淘汰，而改以MATLAB为平台加以重建。
在时间进入20世纪九十年代的时候，MATLAB已经成为国际控制界公认的标准计算软件。
到九十年代初期，在国际上30几个数学类科技应用软件中，MATLAB在数值计算方面独占鳌头。
2.5.2MATLAB组成
MATLAB系统由五个主要部分组成，下面分别加以介绍：
（1）MATLAB语言体系。
MATLAB是高层次的矩阵/数组语言，具有条件控制、函数调用、数据结构、输入输出、面向对象等程序语言特性。
利用它既可以进行小规模编程，完成算法设计和算法实验的基本任务，也可以进行大规模编程，开发复杂的应用程序。
（2）MATLAB工作环境。
这是对MATLAB提高给用户使用的管理功能的总称。
包括管理工作空间中的变量输入输出的方式和方法，以及开发、调试、管理M文件的各种工具。
（3）图形句相系统。
这是MATLAB图形系统的基础，包括完成2D和3D数据图示、图象处理、动画生成、图形显示等功能的高层MATLAB命令，也包括用户对图形图象等对象进行野性控制的低层MATLAB命令，以及开发GUI应用程序的各种工具。
（4）MATLAB数学函数库。
这是对MATLAB使用的各种数学算法
的总称。
包括各种初等函数的算法，也包括矩阵运算、矩阵分析等高层次数学算法。
（5）MATLAB应用程序接口（API）。
这是MATLAB为用户提供的一个函数库，使得用户能够在MATLAB环境中使用C程序或FORTRAN程序，包括从MATLAB中调用程序（动态连接），读写MAT文件的功能。
可以看出MATLAB是一个功能十分强大的系统，是集数值计算、图形管理、程序开发为一体的环境。
除此之外，MATLAB还具有很强的功能扩展能力，与它的系统一起，可以配备各种各样的工具箱。
3频谱分析与算例
3.1声音信号频谱分析
MATLAB语言是一种数据分析和处理功能十分强大的计算机应用软件，它可以将声音文件变换为离散的数据文件，然后利用其强大的矩阵运算能力处理数据，如数字滤波、傅里叶变换、时域和频域分析、声音回放以及各种分析图的呈现等等。
下面以语音信号的波形图、频谱图分析为例来说明MATLAB在语音信号处理中的具体实现方法。
图3-1声音信号分析图
程序代码：
[y0,fs0,nbits0]=wavread（'D:
\work\song.wav'）;
sound（y0,fs0,nbits0）;
n0=length（y0）;
s0=y0;
S0=fft（s0）;
subplot（2,1,1）;
plot（s0）;
title（'信号波形'）;
grid;
subplot（2,1,2）;
plot（abs（S0）,'g'）;
title（'信号频谱'）;
grid;
上程序将语音文件“M.wav”进行频谱分析，分析结果如图3-1所示。
3.2图像信号频谱分析
纹理图像的频谱可以通过离散傅里叶变换（DFT）得到。
用
表示一幅空域纹理图像，用
表示该图像的频谱，图像的大小为M×N，则
和
质检可以通过DFT计算，计算公式如下：
其中能量谱可采用公式：
基于傅立叶能量谱的纹理图像分析的前提是假设纹理有不同的正弦波组成。
理想正弦分布的纹理图像，是最为典型的纹理图像之一，下面讨论理想正弦分布的纹理图像的仿真及其频谱特征分析。
3-2原始图像
利用FFT算法对上面图像信号进行频谱分析。
程序代码为：
I=imread（'1.tif'）
I=rgb2gray（I）;
imshow（I）;
fftI=fft2（I）;
sfftI=fftshift（fftI）;
RR=real（sfftI）;
II=imag（sfftI）;
A=sqrt（RR.^2+II.^2）;
A=（A-min（min（A）））/（max（max（A））-min（min（A）））*225;
figure;
imshow（A）;
在MATLAB中执行了FFT后，使用了fftshift函数调整，以使频谱图像的原点从起始点（0，0），移到图像的中心点（M/2,N/2），对应的图3-3的傅立叶频谱能量图。
从图中可以看出：
竖直方向理想单一频率的正弦分布纹理的频谱能量集中在水平方向的三个点。
图3-3纹理频谱图
4有噪声的语音信号分析与去噪
4.1有噪语音信号提取
4.1.1语音信号的采集
利用PC机上的声卡和WINDOWS操作系统可以进行数字信号的采集。
将话筒输入计算机的语音输入插口上,启动录音机。
按下录音按钮,接着对话筒说话“语音信号处理”,说完后停止录音,屏幕左侧将显示所录声音的长度。
点击放音按钮,可以实现所录音的重现。
以文件名“speech”保存入C:
\MATLAB6p5\work中。
可以看到,文件存储器的后缀默认为.wav,这是WINDOWS操作系统规定的声音文件存的标准。
4.1.2语音信号的时频分析
MATLAB软件平台下，利用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中。
对语音信号speech.wav进行采样其程序如下：
[y,fs,nbits]=wavered（'speech'）;%把语音信号进行加载入MATLAB仿真软件平台中
fs=
44100
nbits=
16
首先画出语音信号的时域波形，然后对语音信号进行频谱分析。
在MATLAB中利用fft对信号进行快速傅里叶变换，得到信号的频谱特性。
其程序如下：
[y,fs,nbits]=wavread（'speech'）;
sound（y,fs,nbits）;%回放语音信号
n=length（y）;%求出语音信号的长度
Y=fft（y,n）;%傅里叶变换
subplot（2,1,1）;plot（y）;title（'原始信号波形'）;
subplot（2,1,2）;plot（abs（Y））;title（'原始信号频谱'）
程序结果如图4-1所示。
图4-1原始信号特征
4.1.3语音信号加噪与频谱分析
利用MATLAB中的随机函数（rand或randn）产生噪声加入到语音信号中，模仿语音信号被污染，并对其频谱分析。
其程序如下：
[y,fs,nbits]=wavread（'speech'）;
n=length（y）;%求出语音信号的长度
noise=0.21*randn（n,2）;%随机函数产生噪声
s=y+noise;%语音信号加入噪声
subplot（2,1,1）;
plot（s）;title（'加噪语音信号的时域波形'）;
S=fft（s）;%傅里叶变换
subplot（2,1,2）;
plot（abs（S））;title（'加噪语音信号的频域波形'）
程序结果如图4-2所示。
图4-2加噪信号特征
4.2设计FIR和IIR数字滤波器
IIR滤波器和FIR滤波器的设计方法完全不同。
IIR滤波器设计方法有间接法和直接法，间接法是借助于模拟滤波器的设计方法进行的。
其设计步骤是：
先设计过渡模拟滤波器得到系统函数H（s），然后将H（s）按某种方法转换成数字滤波器的系统函数H（z）。
FIR滤波器采用间接法，常用的方法有窗函数法、频率采样发和切比雪夫等波纹逼近法。
对于线性相位滤波器，经常采用FIR滤波器。
对于数字带通滤波器的设计，通用方法为双线性变换法。
可以借助于模拟滤波器的频率转换设计一个所需类型的过渡模拟滤波器，再经过双线性变换将其转换策划那个所需的数字滤波器。
具体设计步骤如下：
（1）确定所需类型数字滤波器的技术指标。
（2）将所需类型数字滤波器的边界频率转换成相应的模拟滤波器的边界频率，转换公式为Ω=2/Ttan（0.5ω）
（3）将相应类型的模拟滤波器技术指标转换成模拟低通滤波器技术指标。
（4）设计模拟低通滤波器。
（5）通过频率变换将模拟低通转换成相应类型的过渡模拟滤波器。
（6）采用双线性变换法将相应类型的过渡模拟滤波器转换成所需类型的数字滤波器。
我们知道，脉冲响应不变法的主要缺点是会产生频谱混叠现象，使数字滤波器的频响偏离模拟滤波器的频响特性。
为了克服之一缺点，可以采用双线性变换法。
下面我们总结一下利用模拟滤波器设计IIR数字低通滤波器的步骤：
（1）确定数字低通滤波器的技术指标：
通带边界频率、通带最大衰减，阻带截止频率、阻带最小衰减。
（2）将数字低通滤波器的技术指标转换成相应的模拟低通滤波器的技术指标。
（3）按照模拟低通滤波器的技术指标设计及过渡模拟低通滤波器。
（4）用双线性变换法，模拟滤波器系统函数转换成数字低通滤波器系统函数。
如前所述，IIR滤波器和FIR滤波器的设计方法有很大的区别。
下面我们着重介绍用窗函数法设计FIR滤波器的步骤。
如下：
（1）根据对阻带衰减及过渡带的指标要求，选择串窗数类型（矩形窗、三角窗、汉宁窗、哈明窗、凯塞窗等），并估计窗口长度N。
先按照阻带衰减选择窗函数类型。
原则是在保证阻带衰减满足要求的情况下，尽量选择主瓣的窗函数。
（2）构造希望逼近的频率响应函数。
（3）计算h（n）。
（4）加窗得到设计结果。
接下来，我们根据语音信号的特点给出低通滤波器的性能指标：
fp=1000Hz，fc=1200Hz，As=50db,Ap=1dB
在MATLAB中，可以利用函数fir1设计FIR滤波器，利用函数butter,cheby1和ellip设计IIR滤波器，利用MATLAB中的函数freqz画出各步步器的频率响应。
hn=fir1（M，wc，window），可以指定窗函数向量window。
如果缺省window参数，则fir1默认为哈明窗。
其中可选的窗函数有RectangularBarlrttHammingHannBlackman窗，其相应的都有实现函数。
MATLAB信号处理工具箱函数buttpbuttorbutter是巴特沃斯滤波器设计函数，其有5种调用格式，本课程设计中用到的是[N,wc]=butter（N,wc,Rp,As,’s’）,该格式用于计算巴特沃斯模拟滤波器的阶数N和3dB截止频率wc。
MATLAB信号处理工具箱函数cheblap,cheblord和cheeby1是切比雪夫I型滤波器设计函数。
我们用到的是cheeby1函数，其调用格式如下：
[B,A]=cheby1（N,Rp,wpo,’ftypr’）
[B,A]=cheby1（N,Rp,wpo,’ftypr’,’s’）
函数butter,cheby1和ellip设计IIR滤波器时都是默认的双线性变换法，所以在设计滤波器时只需要代入相应的实现函数即可。
下面我们将给出FIR和IIR数字滤波器的主要程序。
IIR低通滤波器程序：
Ft=8000;
Fp=1000;
Fs=1200;
wp=2*pi*Fp/Ft;
ws=2*pi*Fs/Ft;
fp=2*Ft*tan（wp/2）;
fs=2*Fs*tan（wp/2）;
[n11,wn11]=buttord（wp,ws,1,50,'s'）;%求低通滤波器的阶数和截止频率
[b11,a11]=butter（n11,wn11,'s'）;%求S域的频率响应的参数
[num11,den11]=bilinear（b11,a11,0.5）;%利用双线性变换实现频率响应S域到Z域的变换
[h,w]=freqz（num11,den11）;%根据参数求出频率响应
plot（w*8000*0.5/pi,abs（h））;
legend（'用butter设计'）;
grid
生成相应图片如图4-3所示。
图4-3IIR低通滤波器
相应的IIR带通滤波器程序见附录1
生成相应图片见图4-4所示。
图4-4IIR带通滤波器
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,ma
升级会员 