基于MATLAB的简易声音信号频谱分析仪设计

基于MATLAB的简易声音信号频谱分析仪设计摘要语音信号处理技术是语音处理领域中新近发展起来的一个学科分支，而频谱分析技术是进行语音信号处理的基础。DFT及FFT变换是进行数字信号频谱分析的重要方法。DFT是FFT的基础,FFT是DFT的快速算法。MATLAB是一个数据分析和处理功能十分强大的工程实用软件，运用它来进行语音信号的采集、分析和处理相当方便。本文介绍了在MATLAB环境中如何采集声音信号和采集后的频谱分析方法，并使用MATLAB软件的GUI模块，设计了一个简易的声音信号频谱分析仪。关键字：MATLAB，FFT，声音信号，频谱分析1概述随着信息时代和数字世界的到来,数字信号处理己成为当今一门极其重要的学科和技术领域,数字信号处理在通信、语音、图像、自动控制、医疗和家用电器等众多领域得到了广泛的应用。任意一个信号都具有时域与频域特性,信号的频谱完全代表了信号,因而研究信号的频谱就等于研究信号本身。通常从频域角度对信号进行分析与处理,容易对信号的特性获得深入的了解。因此,信号的频谱分析是数字信号处理技术中的一种较为重要的工具【2】。声卡是计算机最基本的配置硬件之一，价格便宜，使用方便。MATLAB工具箱集成了一些语音处理功能函数。本文将给出基于声卡与MATLAB的声音信号频谱分析仪的设计原理与实现方法。2设计原理频谱分析用傅立叶变换将波形x(t)变换为频谱X(f)，从另一角度来了解信号特征。常见傅里叶变换有DFT和FFT。DFT是FFT的基础,FFT是DFT的快速算法,在MATLAB中可以利用函数fft来计算序列的离散傅里叶变换DFT。FFT是时域和频域转换的基本运算。2.1离散傅里叶级数如果x(n)表示周期为N的周期序列，即：()()xnxnkNk为任意整数（2-1）周期序列用离散的傅里叶级数来表达，其表达式如下：1(2/)01()()NjNknkxnXkeN（2-2）式（2-2）称为周期序列的离散傅里叶变换的级数表示。对上式进行离散傅里叶逆变换，得：1(2/)0()()NjNknnXkxne（2-3）式（2-3）称为周期序列的离散傅里叶逆变换的级数表示。记：(2/)jNNWe（2-4）这样，结合式（2-2）和式（2-3）周期序列的离散傅里叶级数对可表示为：1010()()1()()NknNnNknNkXkxnWxnXkWN（2-5）2.2DFT和FFT变换对于给定的一段时域信号，可以通过傅里叶变换得到相应的频域信号。计算公式如下：1100()()cos(2)()sin(2)NNnnXfxnfntjxnfntt（2-6）上式中，N为样本点数，1/tFs为采样时间间隔。采样信号的频谱是一个连续的频谱，故采用离散傅里叶变换（即DFT），计算公式如下：1(2/)0()(),0,1,2,,1NjNknnXkfxnetkN（2-7）上式中，N为样本点数，1/tFs为采样时间间隔，/fFsN。由于采用式（2-7）进行计算时，有大量的指数（等价于三角函数）运算，效率很低，因此实际中多采用快速傅里叶变换（即FFT）。其原理是通过选择和重新排列，将重复的三角函数计算得到的中间结果保存起来，以减少重复计算带来的时间浪费。由于三角函数计算的重复量相当大，故FFT能极大地提高运算效率。2.3汉宁窗频谱修正采用FFT算法计算信号频谱时，设数据点数为N，采样频率为Fs。则计算得到的离散频率点为：/,0,1,2,,/2ifiFsNiN（2-8）如果信号中的频率分量与频率取样点不重合，则只能按四舍五入的原则，取相邻的频率取样点谱线值代替。这就产生了栅栏效应。频谱的离散取样造成了栅栏效应，谱峰越尖锐，产生误差的可能性就越大。实际应用中，由于信号截断的原因，产生了能量泄漏，即使信号频率与频谱离散取样点不相等，也能得到该频率分量的一个近似值。信号截断带来的能量泄漏分主瓣泄漏和旁瓣泄漏，主瓣泄漏可以减小因栅栏效应带来的谱峰幅值估计误差，而旁瓣泄漏则是完全有害的。实际应用时，可以通过使用截断函数(窗函数)来减小栅栏效应。下面仅以汉宁窗函数为例，说明其工作原理。汉宁窗函数是余弦平方函数，又称为升余弦函数，它的时域形式可以表示为：11cos,2()0,ttTTTwttT（2-9）它的频率幅度特性函数为：sinsinsin1()2TTTWTTT（2-10）汉宁窗的时域和频域曲线图如下所示：汉宁窗时域波形曲线图汉宁窗频域特性曲线图在MATLAB中，生成汉宁窗的函数是hanning。使用该函数进行频谱修正时，先生成一个和待修正的样本具有相同点数的汉宁窗。然后，将原样本序列和生成的汉宁窗序列相乘，得到修正后的样本。最后，对修正后的样本进行FFT变换，即可得到修正后的频谱特性曲线。3MATLAB程序设计3.1图形界面设计首先打开MATLAB，在命令窗口中输入guide命令进入GUI图形设计界面。再新建一个空白的图形界面文件，添加如下的控件并设计它们的布局。（1）添加3个axes控件，用于显示时域波形图和频域频谱图；（2）添加7个statictext控件，用于窗口及其他控件的说明使用；（3）添加4个panel控件，将一组相关的控件框在一起；（4）添加5个edit控件，用于输入和显示幅值、频率等参数值；（5）添加1个pop-upmenu控件，用于选择信号发生器产生的信号类型；（6）添加11个pushbutton控件，其中3个用于控制输出显示相应的信号波形和频谱，2个用于控制播放声音信号，其余6个用于控制3个坐标轴的放大和缩小。双击各个控件，打开其属性编辑窗，即可修改其名称、颜色、大小、初始值及位置等属性。最终编辑好的界面如下图所示：3.2MATLAB编程当界面控件及布局创建完成以后，点击运行即可自动生成包含各控件回调函数在内的m文件。MATLAB对于输入框、按钮及滑动条等控件的响应都是通过自动调用相应的回调函数来实现的。回调函数即在一定的操作下自动执行的指令代码。本次简易声音信号频谱分析仪设计的主要功能有声音文件的打开读取，声音信号的采集录制，信号发生器产生波形，时域信号的频谱分析以及频谱的汉宁窗口校正。下面就依次介绍实现各功能M代码的编写。（1）.wav声音文件的打开读取。实现的代码及说明如下：functionPB_Open_Callback(hObject,eventdata,handles)[filename]=uigetfile('*.wav','选择声音文件');[y,fs]=wavread(filename);handles.y=y;N=size(handles.y);guidata(hObject,handles);t=0:1/fs:(N(1)-1)/fs;plot(handles.axes1,t,handles.y);xlabel(handles.axes1,'Time(s)','fontweight','bold');ylabel(handles.axes1,'Amplitude','fontweight','bold');grid(handles.axes1);上面代码为文件打开按钮的回调函数中的一部分，第一句打开文件对话框，限定选择.wav文件，返回选择的文件名。第二句读取打开的声音文件，并获取音频采样率的值。接着将获得的信号数据存入handles句柄。然后，根据获取到的音频采样率和数据长度还原出时间轴序列。最后将信号波形输出到axes1坐标轴上。（2）声音信号的采集录制。实现的代码如下：functionPB_Record_Callback(hObject,eventdata,handles)fs=str2double(get(handles.edit_Fs,'String'));%获取采样频率值handles.y=wavrecord(5*fs,fs,'int16');%录制声音，并设定时间guidata(hObject,handles);上面代码为录制声音按钮的回调函数其中的一部分。首先利用get函数获取采样率编辑框中的参数值，然后用该采样率录制一段声音信号，并将它保存到handles数据中，方便后面的信号处理和播放。（3）信号发生器产生波形。实现的代码如下：functionPB_SignalGeneration_Callback(hObject,eventdata,handles)fs=str2double(get(handles.edit_Fs,'String'));a=str2double(get(handles.edit_amplitude,'String'));f=str2double(get(handles.edit_frequency,'String'));t=0:1/fs:1.0;i=get(handles.PM_Signaltype,'Value');switchicase1handles.y=a*sin(2*pi*f*t);case2handles.y=a*square(2*pi*f*t);case3handles.y=a*sawtooth(2*pi*f*t,0.5);case4handles.y=a*sawtooth(2*pi*f*t);case5handles.y=sqrt(a)*randn(size(t));endguidata(hObject,handles);上面代码为信号产生分析按钮的回调函数其中的一部分。前三句利用get函数获取到采样频率、产生信号的幅度和频率值，然后再获取弹出菜单控件的选择项序号，根据选择项序号的值，用switch语句和MATLAB自带的波形函数产生相对应的信号波形数据。最后将生成的波形数据保存到handles数据中，方便后面的信号处理。（4）时域信号的频谱分析。实现的代码如下：sample=handles.y(1:N);%设置要进行傅里叶变换的样本sample=double(sample);P=fft(sample,N);%对样本数据进行FFT变换Pyy=2*sqrt(P.*conj(P))/N;%计算频谱幅度f=linspace(0,fs/2,N/2);%还原真实频率plot(handles.axes2,f,Pyy(1:N/2));%画频谱图[Y,I]=max(Pyy(1:N/2));%获取频谱图中的谱峰频率和幅值set(handles.edit_TestA,'string',num2str(Y));set(handles.edit_TestF,'string',num2str(I));xlabel(handles.axes2,'Freqency(Hz)','fontweight','bold');ylabel(handles.axes2,'Amplitude','fontweight','bold');grid(handles.axes2);上面代码为对已存好的信号时域数据进行频域分析的实现过程。首先，截取时域信号中的一段数据作为进行傅里叶变换的样本数据。再用MATLAB自带的fft函数进行FFT变换。接着计算频谱幅度，并获取谱峰的频率和幅度。最后，绘出频谱图并显示谱峰参数。（5）频谱的汉宁窗口校正。实现的代码如下：h=hanning(N);%生成一个汉宁窗口yn=sample.*h';%用汉宁窗修正样本数据Pn=fft(yn,N);%对修正后的样本数据进行FFT变换Pyyn=2*sqrt(Pn.*conj(Pn))/N;plot(handles.axes3,f,Pyyn(1:N/2));xlabel(handles.axes3,'Freqency(Hz)','fontweight','bold');ylabel(handles.axes3,'Amplitude','fontweight','bold');legend(handles.axes3,'