语音信号FFT分析.doc

播码按键语音信号的FFT分析摘要：FFT也叫快速傅氏变换，是离散傅氏变换的快速算法，它是根据离散傅氏变换的奇、偶、虚、实等特性，对离散傅立叶变换的算法进行改进获得。它对傅氏变换的理论并没有新的发现，但是对于在计算机系统或者说数字系统中应用离散傅立叶变换，可以说是进了一大步。本文简要分析了FFT算法的过程以及对电话按键音的不同频率的时域分析和频域分析波形的应用实例。关键词：FFT；MATLAB;按键音；时域；频域；频谱结构1、 FFT的概念FFT也叫快速傅里叶变换，是离散傅里叶变换的快速算法。2、 FFT的工作原理由于严格按照数学推导来阐述FFT的运算原理会带来许多公式编辑的困难，在这里简要的用文字语言对按时间抽取法的基2算法加以说明。DFT的运算中，是利用序列和WNnk按照N为序列索引分别加权求和而得到的。我们单独分析WNnk的 特性可以看出其具有对称性和周期性。因此，我们利用这种性质可以将DFT变换成一系列的迭代运算来完成整个DFT运算过程。整个过程中：1、 我们利用其奇偶性将一个N点DFT运算分成N/2点DFT运算，这样便完成了一个N点DFT的前一半DFT运算。2、 一个N点DFT运算的后一半运算可以利用系数的周期性： r（k+N/2）=WN/2rk和算式 N/2+K=k来推导完成。最后可以得到一个N点DFT的全部运算式： k k其中: k=0,1,2.N/21 3、分解后的序列做DFT运算时任然可以继续按照这种方法再次分为N/2点DFT运算，如果按照这种方法继续进行下去，最终将得到一个运算单元只含2点的DFT运算。 三、 FFT节约的计算量计算 在按时间抽取法中，对于一个N点的序列，所要完成的乘法次数为： 要完成的加法次数为： 而直接DFT运算复乘法数为N2，假设计算1024点DFT时，则FFT和DFT计算的复乘法数目比值是205。因此，FFT算法节约了许多的运算量。4、 一串播码按键语音信号的FFT分析1、按键音的来源在本次实验中，采用的被分析按键音来自于普通手机的播码音，通过MP4录制后保存为WAV的音频格式，然后上传至PC机。由于WAV格式又分为PCM和IMA ADPCM两种不同的格式，MATLAB中wavread函数只能对PCM格式的wav进行处理，而我用的MP4恰好是非PCM格式的wav文件，因此，由MP4录制的语音经过Adensoft Audio MP3 converter格式转换器转换后得到MATLAB能打开的PCM格式WAV文件。将其命名为“q”。2、 语音文件的打开方法 将得到的命名为“q”的WAV文件保存在MATLAB文件下的WORK文件夹以下，然后用函数： y,Fs,bits=wavread(q.wav)能够打开这一串录制好的语音信号。其中，y代表着语音信号的变量名；FS是采样率，BITS是采样位数。q.wav是保存的源语音文件。wavread是函数名。当然，我们也可以回放这一串语音信号，通过调用sound函数即可： sound(y,Fs,bits)3、 语音信号的时域波形显示通过调用MATLAB中的显示命令可以将这串语音信号的时域波形显示出来。通过编辑，从读取语音信号到时域显示部分的MATLAB语言程序如下：y,Fs,bits=wavread(q.wav)；%读取WAV语音文件并命名为y；sound(y,Fs,bits)；%回放语音一次；figure；%产生一个新界面；plot（y）；%绘制时域波形；注：由于这里的输入法不同，所产生的ASICALL码是不一样的，因此复制到MATLAB里面是不被通过的，应重新输入才可。通过运算，得到的原始时域波形如下图所示：此时横坐标的长度是语音信号的名称y的数据序列长度，纵坐标的值是语音信号y的每一个序列值所对应的值得大小。其中可以看出录音一共产生了10次幅度突升的大集合，其实也正对应着0,1,2,3、9等10个数字键的按入。最后的一点幅度是关闭按钮时产生的噪音，可能会对分析的频谱产生噪音污染，但其总体功率不大，任然可以分辨出来。 录音信号的时域波形采用以下程序处理后可以把横轴时间单位改为s： n=length(y); %取y数据的长度； t=(0:n-1)/Fs; %求时间单位s的运算公式； plot(t,y)；xlabel（s）；title（时域波形图）%绘制以时间单位s的时域波形；这样处理后的波形如下图所示： 换算为时间单位的时域波形图4、 语音信号的频域分析在这里，我们采用快速傅里叶变换即FFT分析产生频谱结构的可视化界面，只需要调用相应的函数fft即可： Y=fft（y）；做频谱分析并且显示的程序如下：x=fft（y）；%完成FFT运算，x即运算后的数组；Y=abs（x）；%对FFT运算后的数组x取绝对值，得到幅度值；figure；%产生新界面；plot（Y）；title（原始信号幅度值）；%绘制幅度值图形；经过分析后的时域幅度值波形如下图所示： 录音信号的频域波形图由图中可知，此时得到的频域波形图是对称的并且由于全部显示没有拉开间距，经过以下程序处理后得到的波形将有所改观：W=fftshift（Y）；%移位处理,然后以0为中心；f=linspace（Fs/2，Fs/2，n）；%设置显示范围并且将频域横轴改为以Hz为单位；plot（f，W）；xlabel（频率（Hz）；title（频域波形图）；%绘制频域图形；此时图形如下所示： 单位为Hz的频域波形图由图中可知频率分量主要集中在500Hz和2KHz左右。通过以下显示程序处理后会看到详细的频谱峰值，该显示处理程序主要是调用axis函数：axis（100，2000，0，1600）;%设置横轴最小显示值和最大显示值分别为100Hz，2000Hz；纵轴最小显示值和最大显示值分别为0,1600；截取波形图如下： 缩进频域波形图由图中可以看出，目测峰值频率分量主要集中在700Hz、770Hz、850Hz、950Hz、1200Hz、1330Hz、1470Hz左右，共7个峰值点。这和双音多频数据表上的频率组合几乎完全重合。由于该串语音信号没有录制ABCD等按键，因此这里没有出现1630Hz左右的峰值点。这是通过图中的波形凭观测读取的，是一个大概值。这样便对完整的分析一串语音信号的频谱建立起了程序分析基础。下面，我们分别分析09号键的频谱成分。五、单个按键语音的频谱结构分析1、1号键频谱结构经过FFT运算，然后用语句： axis(650,1300,0,1600)加大缩进量后可以看到1号键的频谱结构如下图所示：由图中可以看到其频谱结构主要由700Hz和1200Hz组成，这和理论值是非常接近的。如果用语句： axis(697,1209,0,1600)的话，将看到与理论值完整的重合。但由于这样的截图不太直观，因此没有截下来，下面也是一样不再阐述。 1号键的频谱结构2、2号键的频谱结构经过FFT运算后然后用语句：axis(680,1350,0,1600)加大缩进量后得到的2号键频谱结构如下图所示： 2号键的频谱结构 由图中可以看到2号键的频谱结构是由695Hz和1340Hz组成，和理论值697Hz和1336Hz也非常的接近。3、3号键的频谱结构经过FFT运算后然后用语句：axis(680,1500,0,1600)加大缩进量后得到的3号键频谱结构如下图所示： 3号键的频谱结构由图中可以看到3号键的频谱结构由700Hz和1480Hz构成，这只是目测数据具有误差的，和理论值697Hz和1477Hz已经非常的接近了。4、4号键的频谱结构经过FFT运算后然后用语句： axis(750,1220,0,1600)加大缩进量后得到的4号键频谱结构如下图所示： 4号键的频谱结构 由图中可以看到4号键的频谱结构由770Hz和1270Hz组成，和理论值770Hz和1209Hz已经很接近。5、5号键频谱结构经过FFT运算后然后用语句： axis(750,1350,0,1600)加大缩进量后得到的5号键频谱结构如下图所示： 5号键频谱结构 由图中可以看到5号键频谱结构主要由770Hz和1340Hz组成，和理论值770Hz和1336Hz很接近。6、6号键的频谱结构经过FFT运算后然后用语句： axis(750,1500,0,1500)加大缩进量后得到的6号键频谱结构如下图所示： 6号键的频谱结构由图中可以看到6号键的频谱结构主要由770Hz和1480Hz构成，和理论值770Hz和1477Hz也非常的接近。7、7号键的频谱结构经过FFT运算后然后用语句： axis(800,1250,0,1000)加大缩进量后得到的7号键频谱结构如下图所示： 7号键的频谱结构由图中可以看到7号键的频谱结构由855Hz和1210Hz组成，这和理论值852Hz和1209Hz很接近。8、8号键的频谱结构经过FFT运算后然后用语句：axis(800,1350,0,1600)加大缩进量后得到的8号键频谱结构如下图所示： 8号键的频谱结构由图中可以看到 8号键的频谱结构主要由855Hz和1340Hz组成，和理论值852和1336Hz很接近。9、9号键的频谱结构经过FFT运算后然后用语句：axis(800,1500,0,1100)加大缩进量后得到的9号键频谱结构如下图所示： 9号键的频谱结构由图中可以看到9号键的频谱结构主要由850Hz和1480Hz组成，和理论值852Hz和1477Hz很接近。10、0号键的频谱结构经过FFT运算后然后用语句：axis(900,1350,0,600)加大缩进量后得到的0号键频谱结构如下图所示： 0号键的频谱结构由图中可以看到0号键的频谱结构主要是940Hz和1340Hz组成的，和理论值941Hz和1336Hz很接近。以上便是全部的双音多频从09号按键的频谱分析。 六、FFT分析程序 这里将所有的分析程序附在下面，其可以分析一串或者单个的语音信号。单个双音多频按键的缩进语句在上面已经阐述，以下程序都采用统一的视野缩进显示。另外，以下程序输入法是经过处理的，可以通过直接复制到MATLAB界面运行。如果要分析具体的一段语音，可以把语音文件名改为“q”，或者把程序第一条wavread语句里面的“q”改为相应的语音文件名即可：y,Fs,bits=wavread(q.wav);sound(y,Fs,bits);figure;plot(y);n=length(y);t=(0:n-1)/Fs;figure;plot(t,y);xlabel(s);title(时域波形);x=fft(y);Y=abs(x);figure;plot(Y);title(原始信号幅度值);W=fftshift(Y);figure;plot(W);title(中心频率点移位处理后的频谱图形);f=linspace(-Fs/2