信号处理实验七音频频谱分析仪设计及实现

1、哈尔滨工程大学实 验 报 告实 验 名 称： 离散时间滤波器设计 班 级： 电子信息工程4班 学 号： 姓 名： 实 验 时 间： 2016年10月31日18：30 成 绩：_指 导 教 师： 栾晓明 实验室名称： 数字信号处理实验室 哈尔滨工程大学实验室与资产管理处 制实验七 音频频谱分析仪设计与实现一、 实验原理MATLAB是一个数据分析和处理功能十分强大的工程实用软件，其数据采集工具箱为实现数据的输入和输出提供了十分方便的函数命令。本实验要求基于声卡和MTLAB实现音频信号频谱分析仪的设计原理与实现，功能包括：(1)音频信号输入，从声卡输入、从WAV文件输入、从标准信号发生器输入；(2)

2、信号波形分析，包括幅值、频率、周期、相位的估计、以及统计量峰值、均值、均方值和方差的计算。(3)信号频谱分析，频率、周期的统计，同行显示幅值谱、相位谱、实频谱、虚频谱和功率谱的曲线。1、频率(周期)检测对周期信号来说，可以用时域波形分析来确定信号的周期，也就是计算相邻的两个信号波峰的时间差、或过零点的时间差。这里采用过零点(ti)的时间差T(周期)。频率即为f = 1/T，由于能够求得多个T值(ti有多个)，故采用它们的平均值作为周期的估计值。2、幅值检测在一个周期内，求出信号最大值ymax与最小值ymin的差的一半，即A = (ymax - ymin)/2，同样，也会求出多个A值，但第1个A

3、值对应的ymax和ymin不是在一个周期内搜索得到的，故以除第1个以外的A值的平均作为幅值的估计值。3、相位检测采用过零法，即通过判断与同频零相位信号过零点时刻，计算其时间差，然后换成相应的相位差。=2(1-ti/T)，x表示x的小数部分，同样，以的平均值作为相位的估计值。频率、幅值和相位估计的流程如图1所示。4、数字信号统计量估计(1) 峰值P的估计在样本数据x中找出最大值与最小值，其差值为双峰值，双峰值的一半即为峰值。P=0.5max(yi)-min(yi)(2)均值估计 式中，N为样本容量，下同。(3) 均方值估计 (4) 方差估计 开始n=1,ymax =ymin=0i=3Yyi=10

4、且yi=2=0且yi+10NNymaxN-1?i+NYT为（ti -ti-1)的平均值，i=1.2,f=1/Tn+ymin =yiA 为Ai的平均值，i=2,3,yminyi为2（1-ti /Ti)的平均值结束图 1频率,幅值和相位估计的流程图其中tin表示第n个过零点，yi为第i个采样点的值，Fs为采样频率。5、频谱分析原理时域分析只能反映信号的幅值随时间的变化情况，除单频率分量的简单波形外，很难明确提示信号的频率组成和各频率分量大小，而频谱分析能很好的解决此问题。由于从频域能获得的主要是频率信息，所以本节主要介绍频率(周期)的估计与频谱图的生成。(1)DFT与FFT对于给定的时域信号y，可

5、以通过Fourier变换得到频域信息Y，Y可按下式计算式中，N为样本容量，t = 1/Fs为采样间隔。采样信号的频谱是一个连续的频谱，不可能计算出所有的点的值，故采用离散Fourier变换(DFT)，即式中，f = Fs/N。但上式的计算效率很低，因为有大量的指数(等价于三角函数)运算，故实际中多采用快速Fourier变换(FFT)。其原理即是将重复的三角函数算计的中间结果保存起来，以减少重复三角函数计算带来的时间浪费。由于三角函数计算的重复量相当大，故FFT能极大地提高运算效率。(2)频率、周期的估计对于Y(kf)，如果当kf = f时，Y(kf)取最大值，则f为频率的估计值，由于采样间隔的

6、误差，f也存在误差，其误差最大为f / 2，周期T=1/f。 从原理上可以看出，如果在标准信号中混有噪声，用上述方法仍能够精确地估计出原标准信号的频率和周期。(3)频谱图为了直观地表示信号的频率特性，工程上常常将Fourier变换的结果用图形的方式表示，即频谱图。以频率f为横坐标，|Y(f)|为纵坐标，可以得到幅值谱；以频率f为横坐标，argY(f)为纵坐标，可以得到相位谱；以频率f为横坐标，ReY(f)为纵坐标，可以得到实频谱；以频率f为横坐标，ImY(f)为纵坐标，可以得到虚频谱。 根据采样定理，只有频率不超过Fs/2的信号才能被正确采集，即Fourier变换的结果中频率大于Fs/2的部分

7、是不正确的部分，故不在频谱图中显示。即横坐标f 0, Fs/2(4)频谱图模块化就是把程序划分成独立命名且可独立访问的模块，每个模块完成一个子功能，把这些模块集成起来构成一个整体，可以完成指定的功能满足用户需求。根据人类解决一般问题的经验，如果一个问题由两个问题组合而成，那么它的复杂程度大于分别考虑每个问题时的复杂程度之和，也就是说把复杂的问题分解成许多容易解决的小问题，原来的问题也就容易解决了。这就是模块化的根据。在模块划分时应遵循如下规则：1.改进软件结构提高模块独立性；2.模块规模应该适中；3.深度、宽度、扇出和扇入都应适当；4.模块的作用域应该在控制域之内；5.力争降低模块接口的复杂程

8、度；6.设计单入口单出口的模块；7.模块功能应该可以预测。本着上述的启发式规则，对软件进行如图 2所示的模块划分。图 2频谱分析仪的模块划分二、 界面设计MATLAB是Mathworks公司推出的数学软件，它将数值分析、矩阵计算、信号处理和图形显示结合在一起，为众多学科领域提供了一种简洁、高效的编程工具。它提供的GUIDE工具为可视化编程工具，使得软件的界面设计像VB一样方便。故本文采用MATLAB作为编程语言实现声音信号频谱分析仪，以下所讲的都是在MATLAB2014a环境中。为了实现预期的功能，设计如图 3所示的界面。 图 3声音信号频谱分析仪最上面的部分为标题区，用于显示软件标题等信息，

9、不具人机交互功能。标题区下面是信号输入区，包含3种输入方式：声卡输入，WAV文件输入，信号发生器输入。考虑到WAV文件可能是多声道，故提供了声道选择的界面，因为每次只能对单个声道进行分析。在信号发生器中加入了混迭选项，从而可以将产生的信号与原有的信号进行混迭。输入方式界面应该具有：只有当每个单选框被选中时才允许使用对应的输入框、按钮等；其中采样点数输入框在声卡与WAV文件的输入方式下作为输出，在信号发生器的输入方式下作为输入。输入方式单选框程序代码为:(1)声卡单选框程序代码set(findobj(Tag,recordtime),enable,on);h=findobj(Tag,filenam

10、e);set(h,enable,off);h=findobj(Tag,freq);set(h,enable,off);h=findobj(Tag,amp);set(h,enable,off);h=findobj(Tag,phase);set(h,enable,off);set(handles.channel,enable,off);set(handles.fileopen,enable,off);set(handles.gensig,enable,off);set(handles.wavetype,enable,off);set(handles.add,enable,off);set(hand

11、les.startrecord,enable,on);（2）WAV文件单选框程序代码h=findobj(Tag,filename);set(h,enable,on);h=findobj(Tag,freq);set(h,enable,off);h=findobj(Tag,amp);set(h,enable,off);h=findobj(Tag,phase);set(h,enable,off);set(findobj(Tag,recordtime),enable,off);set(handles.channel,enable,on);set(handles.fileopen,enable,on);

12、set(handles.gensig,enable,off);set(handles.wavetype,enable,off);set(handles.add,enable,off);set(handles.startrecord,enable,off);（3）信号发生器单选框程序代码h=findobj(Tag,filename);set(h,enable,off);h=findobj(Tag,freq);set(h,enable,on);h=findobj(Tag,amp);set(h,enable,on);h=findobj(Tag,phase);set(h,enable,on);set(

13、findobj(Tag,recordtime),enable,off);set(handles.channel,enable,off);set(handles.fileopen,enable,off);set(handles.gensig,enable,on);set(handles.wavetype,enable,on);set(handles.add,enable,on);set(handles.startrecord,enable,off); 再往下是分析区，对于WAV文件及录音的信号，有时只对其中一部分信号进行分析，故提供了分析对象范围设定的界面。另外就是时域分析与频域分析的按钮，该软

14、件的核心代码都在这两个按钮的回调函数中。分析区下面是分析结果区，用于显示波形基本参数与统计量的计算结果。分析结果区的下面是波形显示区，用于显示时域波形，在录音结束、打开WAV文件成功或者信号发生器生成波形时会更新显示。右边为频谱图显示区，用于显示各种频谱的谱线，在点击频域分析后会更新显示。1、输入模块的实现采样频率Fs与采样点数N是声音信号输入时共同需要作用的参数，故将其独立出来。下面分别介绍三种输入方式的实现。（1）声卡输入这里声卡输入是指由麦克风录音得到的声音信号的输入，MATLAB提供了wavrecord函数，该函数能够实现读取麦克风录音信号。以下是“开始录音”按钮的回调函数的程序代码。

15、Fs=str2double(get(findobj(Tag,samplerate),String);handles.y=wavrecord(str2double(get(findobj(Tag,recordtime),String)*Fs, Fs,int16);handles.inputtype=1;guidata(hObject,handles);plot(handles.time,handles.y);title(WAVE);ysize=size(handles.y)set(handles.samplenum,String,num2str(ysize(1);（2）WAV文件输入MATLAB

16、提供了wavread函数，该函数能够方便的打开并读取WAV文件中的声音信息，并且同时读取所有声道。以下是“打开文件”按钮回调函数的程序代码。temp = wavread(get(findobj(Tag,filename),String);channel=str2double(get(handles.channel,String);handles.y=temp(:,channel);handles.inputtype=2;guidata(hObject,handles);plot(handles.time,handles.y);title(WAVE);ysize=size(handles.y)s

17、et(handles.samplenum,String,num2str(ysize(1);（3）信号发生器MATLAB有产生标准信号的函数，如sawtooth能够产生三角波或钜齿波，利用get函数获得波形soundtype，频率frequency，幅值amp和相位phase，以下是“生成波形”按钮回调函数程序代码。Fs=str2double(get(findobj(Tag,samplerate),String);N=str2double(get(findobj(Tag,samplenum),String);x=linspace(0,N/Fs,N);soundtype=get(handles.w

18、avetype,Value);frequency=str2double(get(handles.freq,String);amp=str2double(get(handles.amp,String);phase=str2double(get(handles.phase,String);switch soundtypecase 1 y=amp*sin(2*pi*x*frequency+phase);case 2 y=amp*sign(sin(2*pi*x*frequency+phase);case 3 y=amp*sawtooth(2*pi*x*frequency+phase,0.5);case

19、 4 y=amp*sawtooth(2*pi*x*frequency+phase); case 5 y=amp*(2*rand(size(x)-1); otherwise errordlg(Illegal wave type,Choose errer);endif get(handles.add,Value)=0.0 handles.y=y;else handles.y=handles.y+y;endhandles.inputtype=3;guidata(hObject,handles);plot(handles.time,handles.y);title(WAVE);axis(0 N -st

20、r2double(get(handles.amp,String) str2double(get(handles.amp,String);2、分析模块与图形显示模块由于MATLAB的绘图功能很强大，所以图形显示模块不用单独开发，可直接调用plot、axis等函数实现图形显示功能，故图形显示也将在分析模块中给出。（1）时域分析MATLAB提供了mean，std函数，能够方便地计算均值、标准差。以下是“时域分析”按钮回调函数程序代码，其中T为过零检测得到的周期(向量)，amp为过零检测得到的幅值(向量)，n为过零点数。Fs=str2double(get(findobj(Tag,samplerate)

21、,String);N=str2double(get(findobj(Tag,samplenum),String);if handles.inputtype=0 msgbox(No wave exist! Please choose a input type!); return;endn=1;ymax=max(handles.y(1) handles.y(2);ymin=min(handles.y(1) handles.y(2);from=str2double(get(handles.pointfrom,String);to=str2double(get(handles.pointto,Stri

22、ng);if from1 | to-from5; msgbox(Error range!); return;endfor i=from+2:to-1; if handles.y(i-1)0 & handles.y(i-2)=0 & handles.y(i+1)0 if handles.y(i)=0 ti(n)=i; else ti(n)=i-handles.y(i)/(handles.y(i)-handles.y(i-1); end amp(n)=(ymax-ymin)/2; ymax=0; ymin=0; n=n+1; else if ymaxhandles.y(i) ymin=handle

23、s.y(i); end endendn=n-1;for i=1:n-1 T(i)=ti(i+1)-ti(i);endfreq=Fs/mean(T);set(handles.outt,String,1/freq);set(handles.outfreq,String,num2str(freq);set(handles.outamp,String,num2str(mean(amp(2:n-1);phase=2*pi*(1-(ti(1:n-1)-1)./T+floor(ti(1:n-1)-1)./T);set(handles.outphase,String,num2str(mean(phase);s

24、et(handles.outpeak,String,(max(handles.y(from:to)-min(handles.y(from:to)/2);set(handles.outmean,String,mean(handles.y(from:to);set(handles.outmeansquare,String,mean(handles.y(from:to).2);set(handles.outs,String,std(handles.y(from:to)2);（2）频域分析频域分析需要作Fourier变换，MATLAB提供了fft函数，能够方便地实现快速Fourier变换算法。以下是“

25、频域分析”按钮回调函数程序代码。Fs=str2double(get(findobj(Tag,samplerate),String);N=str2double(get(findobj(Tag,samplenum),String);if handles.inputtype=0 msgbox(No wave exist! Please choose a input type!); return;endfrom=str2double(get(handles.pointfrom,String);to=str2double(get(handles.pointto,String);sample=handle

26、s.y(from:to);f=linspace(0,Fs/2,(to-from+1)/2);Y=fft(sample,to-from+1);C,I=max(abs(Y);set(handles.foutt,String,1/f(I);set(handles.foutfreq,String,f(I);Y=Y(1:(to-from+1)/2);plot(handles.plot1,f,2*sqrt(Y.*conj(Y);plot(handles.plot2,f,angle(Y);plot(handles.plot3,f,real(Y);plot(handles.plot4,f,imag(Y);pl

27、ot(handles.plot5,f,abs(Y).2);xlabel(handles.plot1,freqency(Hz);xlabel(handles.plot2,freqency(Hz);xlabel(handles.plot3,freqency(Hz);xlabel(handles.plot4,freqency(Hz);xlabel(handles.plot5,freqency(Hz);ylabel(handles.plot1,amplitude);ylabel(handles.plot2,phase(rad);ylabel(handles.plot3,real);ylabel(han

28、dles.plot4,Imaginary);ylabel(handles.plot5,power);三、 软件运行及结果分析为了分析软件的性能并比较时域分析与频域分析各自的优势，本章给出了两种分析频率估计方法的比较，分析软件的性能在时域和频域的计算精度问题。(1)标准正弦信号的频率估计用信号发生器生成标准正弦信号频率设为100Hz，然后分别进行时域分析与频域分析，得到的结果如图4所示。从图中可以看出，时域分析的结果为f = 100.1965Hz，频域分析的结果为f = 104.49Hz，而标准信号的频率为100Hz，从而对于标准信号时域分析的精度远高于频域分析的精度。(2)带噪声的正弦信号的频率估计先生成幅值100的标准正弦信号，再将幅值100的白噪声信号与其混迭，对最终得到的信号进行时域分析与频域分析，结果如图 5所示，可以看出，时域分析的结果为f