Y
yi=0?
Y
ymax=yi
N
tin=i/Fs
N
Y
tin=[i-yi/（yi-yi-1）]/Fs
An=（ymax-ymin）/2
ymax=ymin=0
i>N-1?
i++
N
Y
T为（ti-ti-1）的平均值，i=1.2,…
f=1/T
n++
ymin=yi
A为Ai的平均值，i=2,3,…
ymin>yi
φ为2π（1-{ti/Ti}）的平均值
结束
图1 频率,幅值和相位估计的流程图
其中tin表示第n个过零点，yi为第i个采样点的值，Fs为采样频率。
5、频谱分析原理
时域分析只能反映信号的幅值随时间的变化情况，除单频率分量的简单波形外，很难明确提示信号的频率组成和各频率分量大小，而频谱分析能很好的解决此问题。
由于从频域能获得的主要是频率信息，所以本节主要介绍频率（周期）的估计与频谱图的生成。
（1）DFT与FFT
对于给定的时域信号y，可以通过Fourier变换得到频域信息Y，Y可按下式计算
式中，N为样本容量，Δt=1/Fs为采样间隔。
采样信号的频谱是一个连续的频谱，不可能计算出所有的点的值，故采用离散Fourier变换（DFT），即
式中，Δf=Fs/N。
但上式的计算效率很低，因为有大量的指数（等价于三角函数）运算，故实际中多采用快速Fourier变换（FFT）。
其原理即是将重复的三角函数算计的中间结果保存起来，以减少重复三角函数计算带来的时间浪费。
由于三角函数计算的重复量相当大，故FFT能极大地提高运算效率。
（2）频率、周期的估计
对于Y（kΔf），如果当kΔf=f时，Y（kΔf）取最大值，则f为频率的估计值，由于采样间隔的误差，f也存在误差，其误差最大为Δf/2，周期T=1/f。
从原理上可以看出，如果在标准信号中混有噪声，用上述方法仍能够精确地估计出原标准信号的频率和周期。
（3）频谱图
为了直观地表示信号的频率特性，工程上常常将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]
（4）频谱图
模块化就是把程序划分成独立命名且可独立访问的模块，每个模块完成一个子功能，把这些模块集成起来构成一个整体，可以完成指定的功能满足用户需求。
根据人类解决一般问题的经验，如果一个问题由两个问题组合而成，那么它的复杂程度大于分别考虑每个问题时的复杂程度之和，也就是说把复杂的问题分解成许多容易解决的小问题，原来的问题也就容易解决了。
这就是模块化的根据。
在模块划分时应遵循如下规则：
1.改进软件结构提高模块独立性；2.模块规模应该适中；3.深度、宽度、扇出和扇入都应适当；4.模块的作用域应该在控制域之内；5.力争降低模块接口的复杂程度；6.设计单入口单出口的模块；7.模块功能应该可以预测。
本着上述的启发式规则，对软件进行如图2所示的模块划分。
图2频谱分析仪的模块划分
二、界面设计
MATLAB是Mathworks公司推出的数学软件，它将数值分析、矩阵计算、信号处理和图形显示结合在一起，为众多学科领域提供了一种简洁、高效的编程工具。
它提供的GUIDE工具为可视化编程工具，使得软件的界面设计像VB一样方便。
故本文采用MATLAB作为编程语言实现声音信号频谱分析仪，以下所讲的都是在MATLAB2014a环境中。
为了实现预期的功能，设计如图3所示的界面。
图3声音信号频谱分析仪
最上面的部分为标题区，用于显示软件标题等信息，不具人机交互功能。
标题区下面是信号输入区，包含3种输入方式：
声卡输入，WAV文件输入，信号发生器输入。
考虑到WAV文件可能是多声道，故提供了声道选择的界面，因为每次只能对单个声道进行分析。
在信号发生器中加入了混迭选项，从而可以将产生的信号与原有的信号进行混迭。
输入方式界面应该具有：
只有当每个单选框被选中时才允许使用对应的输入框、按钮等；其中采样点数输入框在声卡与WAV文件的输入方式下作为输出，在信号发生器的输入方式下作为输入。
输入方式单选框程序代码为:
（1）声卡单选框程序代码
set（findobj（'Tag','recordtime'）,'enable','on'）;
h=findobj（'Tag','filename'）;
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（handles.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'）;
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（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文件及录音的信号，有时只对其中一部分信号进行分析，故提供了分析对象范围设定的界面。
另外就是时域分析与频域分析的按钮，该软件的核心代码都在这两个按钮的回调函数中。
分析区下面是分析结果区，用于显示波形基本参数与统计量的计算结果。
分析结果区的下面是波形显示区，用于显示时域波形，在录音结束、打开WAV文件成功或者信号发生器生成波形时会更新显示。
右边为频谱图显示区，用于显示各种频谱的谱线，在点击频域分析后会更新显示。
1、输入模块的实现
采样频率Fs与采样点数N是声音信号输入时共同需要作用的参数，故将其独立出来。
下面分别介绍三种输入方式的实现。
（1）声卡输入
这里声卡输入是指由麦克风录音得到的声音信号的输入，MATLAB提供了wavrecord函数，该函数能够实现读取麦克风录音信号。
以下是“开始录音”按钮的回调函数的程序代码。
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提供了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）
set（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.wavetype,'Value'）;
frequency=str2double（get（handles.freq,'String'））;
amp=str2double（get（handles.amp,'String'））;
phase=str2double（get（handles.phase,'String'））;
switchsoundtype
case1
y=amp*sin（2*pi*x*frequency+phase）;
case2
y=amp*sign（sin（2*pi*x*frequency+phase））;
case3
y=amp*sawtooth（2*pi*x*frequency+phase,0.5）;
case4
y=amp*sawtooth（2*pi*x*frequency+phase）;
case5
y=amp*（2*rand（size（x））-1）;
otherwise
errordlg（'Illegalwavetype','Chooseerrer'）;
end
ifget（handles.add,'Value'）==0.0
handles.y=y;
else
handles.y=handles.y+y;
end
handles.inputtype=3;
guidata（hObject,handles）;
plot（handles.time,handles.y）;
title（'WAVE'）;
axis（[0N-str2double（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'）,'String'））;
N=str2double（get（findobj（'Tag','samplenum'）,'String'））;
ifhandles.inputtype==0
msgbox（'Nowaveexist!
Pleasechooseainputtype!
'）;
return;
end
n=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,'String'））;
iffrom<1|to-from<5;
msgbox（'Errorrange!
'）;
return;
end
fori=from+2:
to-1;
ifhandles.y（i-1）<0&handles.y（i-2）<0&handles.y（i）>=0&handles.y（i+1）>0
ifhandles.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
ifymaxymax=handles.y（i）;
end
ifymin>handles.y（i）
ymin=handles.y（i）;
end
end
end
n=n-1;
fori=1:
n-1
T（i）=ti（i+1）-ti（i）;
end
freq=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）））;
set（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变换算法。
以下是“频域分析”按钮回调函数程序代码。
Fs=str2double（get（findobj（'Tag','samplerate'）,'String'））;
N=str2double（get（findobj（'Tag','samplenum'）,'String'））;
ifhandles.inputtype==0
msgbox（'Nowaveexist!
Pleasechooseainputtype!
'）;
return;
end
from=str2double（get（handles.pointfrom,'String'））;
to=str2double（get（handles.pointto,'String'））;
sample=handles.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））;
plot（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（handles.plot4,'Imaginary'）;
ylabel（handles.plot5,'power'）;
三、软件运行及结果分析
为了分析软件的性能并比较时域分析与频域分析各自的优势，本章给出了两种分析频率估计方法的比较，分析软件的性能在时域和频域的计算精度问题。
（1）标准正弦信号的频率估计
用信号发生器生成标准正弦信号频率设为100Hz，然后分别进行时域分析与频域分析，得到的结果如图4所示。
从图中可以看出，时域分析的结果为f=100.1965Hz，频域分析的结果为f=104.49Hz，而标准信号的频率为100Hz，从而对于标准信号时域分析的精度远高于频域分析的精度。
（2）带噪声的正弦信号的频率估计
先生成幅值100的标准正弦信号，再将幅值100的白噪声信号与其混迭，对最终得到的信号进行时域分析与频域分析，结果如图5所示，可以看出，时域分析的结果为f=133.824Hz，频域分析的结果为f=100.787Hz，而标准信号的频率为100Hz，从而对于带噪声的正弦信号频域分析的精度远高于时域分析的精度。
结果分析：
在时域，频率估计是使用过零检测的方式计算出，从而对于带噪声的信号既容易造成“误判”，也容易造成“漏判”，且噪声信号越明显，“误判”与“漏判”的可能性越大。
但在没有噪声或噪声很小时，时域分析对每个周期长度的检测是没有累积误差的，故随着样本容量的增大，估计的精度大大提高。
在频域，频率估计是通过找出幅值谱峰值点对
升级会员 