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

1．概述

随着软硬件技术的发展，仪器的智能化与虚拟化已成为未来实验室及研究机构的发展方向^[1]。虚拟仪器技术的优势在于可由用户定义自己的专用仪器系统，且功能灵活，很容易构建，所以应用面极为广泛。基于计算机软硬件平台的虚拟仪器可代替传统的测量仪器，如示波器、逻辑分析仪、信号发生器、频谱分析仪等^[2]。从发展史看，电子测量仪器经历了由模拟仪器、智能仪器到虚拟仪器，由于计算机性能的飞速发展，已把传统仪器远远抛到后面，并给虚拟仪器生产厂家不断带来连锅端的技术更新速率。目前已经有许多较成熟的频谱分析软件，如SpectraLAB、RSAVu、dBFA等。

声卡是多媒体计算机最基本的配置硬件之一，价格便宜，使用方便。MATLAB是一个数据分析和处理功能十分强大的工程实用软件，他的数据采集工具箱为实现数据的输入和输出提供了十分方便的函数和命令^[3]。本文将给出基于声卡与MATLAB的声音信号频谱分析仪的设计原理与实现方法，功能包括：

(1) 音频信号信号输入，从声卡输入、从WAV文件输入、从标准信号发生器输入；

(2) 信号波形分析，包括幅值、频率、周期、相位的估计，以及统计量峰值、均值、均方值和方差的计算；

(3) 信号频谱分析，频率、周期的估计，图形显示幅值谱、相位谱、实频谱、虚频谱和功率谱的曲线。

2．设计原理

2.1波形分析原理

2.1.1 信号频率、幅值和相位估计

(1)频率(周期)检测

对周期信号来说，可以用时域波形分析来确定信号的周期，也就是计算相邻的两个信号波峰的时间差、或过零点的时间差。这里采用过零点(ti)的时间差T(周期)。频率即为f = 1/T，由于能够求得多个T值(ti有多个)，故采用它们的平均值作为周期的估计值。

(2)幅值检测

在一个周期内，求出信号最大值y_max与最小值y_min的差的一半，即A = (y_max - y_min)/2，同样，也会求出多个A值，但第1个A值对应的y_max和y_min不是在一个周期内搜索得到的，故以除第1个以外的A值的平均作为幅值的估计值。

(3)相位检测

采用过零法，即通过判断与同频零相位信号过零点时刻，计算其时间差，然后换成相应的相位差。φ=2π(1-ti/T)，{x}表示x的小数部分，同样，以φ的平均值作为相位的估计值。

频率、幅值和相位估计的流程如图 1所示。

图 1　频率、幅值和相位估计的流程图

其中ti_n表示第n个过零点，y_i为第i个采样点的值，Fs为采样频率。

2.1.2 数字信号统计量估计

(1) 峰值P的估计

在样本数据x中找出最大值与最小值，其差值为双峰值，双峰值的一半即为峰值。

P=0.5[max(y_i)-min(y_i)]

(2)均值估计

式中，N为样本容量，下同。

(3) 均方值估计

(4) 方差估计

 

2.2频谱分析原理

时域分析只能反映信号的幅值随时间的变化情况，除单频率分量的简单波形外，很难明确提示信号的频率组成和各频率分量大小，而频谱分析能很好的解决此问题。由于从频域能获得的主要是频率信息，所以本节主要介绍频率(周期)的估计与频谱图的生成。

2.2.1 DFT与FFT

对于给定的时域信号y，可以通过Fourier变换得到频域信息Y。Y可按下式计算

式中，N为样本容量，Δt = 1/Fs为采样间隔。

采样信号的频谱是一个连续的频谱，不可能计算出所有的点的值，故采用离散Fourier变换(DFT)，即

式中，Δf = Fs/N。但上式的计算效率很低，因为有大量的指数(等价于三角函数)运算，故实际中多采用快速Fourier变换(FFT)。其原理即是将重复的三角函数算计的中间结果保存起来，以减少重复三角函数计算带来的时间浪费。由于三角函数计算的重复量相当大，故FFT能极大地提高运算效率。

2.2.2 频率、周期的估计

对于Y(kΔf)，如果当kΔf = 时，Y(kΔf)取最大值，则为频率的估计值，由于采样间隔的误差，也存在误差，其误差最大为Δf / 2。

周期T=1/f。

从原理上可以看出，如果在标准信号中混有噪声，用上述方法仍能够精确地估计出原标准信号的频率和周期，这个将在下一章做出验证

2.2.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]

2.3. 模块划分

模块化就是把程序划分成独立命名且可独立访问的模块，每个模块完成一个子功能，把这些模块集成起来构成一个整体，可以完成指定的功能满足用户需求。根据人类解决一般问题的经验，如果一个问题由两个问题组合而成，那么它的复杂程度大于分别考虑每个问题时的复杂程度之和，也就是说把复杂的问题分解成许多容易解决的小问题，原来的问题也就容易解决了。这就是模块化的根据。

在模块划分时应遵循如下规则^[4]：改进软件结构提高模块独立性；模块规模应该适中；深度、宽度、扇出和扇入都应适当；模块的作用域应该在控制域之内；力争降低模块接口的复杂程度；设计单入口单出口的模块；模块功能应该可以预测。

本着上述的启发式规则，对软件进行如图 2所示的模块划分。

 

图 2频谱分析仪的模块划分

3．软件实现

3.1界面设计

MATLAB是Mathworks公司推出的数学软件，它将数值分析、矩阵计算、信号处理和图形显示结合在一起，为众多学科领域提供了一种简洁、高效的编程工具。它提供的GUIDE工具为可视化编程工具，使得软件的界面设计像VB一样方便。故本文采用MATLAB作为编程语言实现声音信号频谱分析仪，以下所讲的都是在MATLAB7.0环境中。

为了实现预期的功能，设计如图 3所示的界面。

 

图 3　频谱分析仪的界面设计

最上面的部分为标题区，用于显示软件标题等信息，不具人机交互功能。

再往下是信号输入区，包含3种输入方式，考虑到WAV文件可能是多声道，故提供了声道选择的界面，因为每次只能对单个声道进行分析。在信号发生器中加入了混迭选项，从而可以将产生的信号与原有的信号进行混迭。界面应该具有：只有当每个单选框被选中时才允许使用对应的输入框、按钮等；采样点数输入框在声卡与WAV文件的输入方式下作为输出，在信号发生器的输入方式下作为输入。

再往下是分析区。对于WAV文件及录音的信号，有时只对其中一部分信号进行分析，故提供了分析对象范围设定的界面。另外就是时域分析与频域分析的按钮，该软件的核心代码都在这两个按钮的回调函数中。

分析区下面是分析结果区，用于显示波形基本参数与统计量的计算结果。

分析结果区的下面是波形显示区，用于显示时域波形，在录音结束、打开WAV文件成功或者信号发生器生成波形时会更新显示。

右边为频谱图显示区，用于显示各种频谱的谱线，在点击频域分析后会更新显示。

3.2输入模块的实现

采样频率Fs与采样点数N是声音信号输入时共同需要作用的参数，故将其独立出来。下面为别介绍三种输入方式的实现。

3.2.1 声卡输入

这里声卡输入是指由麦克风录音得到的声音信号的输入，MATLAB提供了wavrecord函数，该函数能够实现读取麦克风录音信号。以下是“开始录音”按钮的回调函数内容。

%首先获得设定的Fs值

Fs=str2double(get(findobj('Tag','samplerate'),'String'));

%根据设定的录音时长进行录音，将其存入handles.y中

handles.y=wavrecord(str2double(get(handles.recordtime,'String'))*Fs, Fs,'int16');

%保存handles结构体，使得handles.y在别的函数中也能使用

guidata(hObject,handles);

%在波形显示区绘出波形

plot(handles.time,handles.y);

title('WAVE');

%将所采到的点的数量输出在“采样点数”中

ysize=size(handles.y)

set(handles.samplenum,'String',num2str(ysize(1)));

3.2.2 WAV文件输入

MATLAB提供了wavread函数，该函数能够方便的打开并读取WAV文件中的声音信息，并且同时读取所有声道。下面是“打开文件”按钮回调函数的部分代码。其它代码与声卡输入的类似。

%从WAV文件中读取的声音信息并临时存放到temp变量中

temp = wavread(get(findobj('Tag','filename'),'String'));

%获得所选择的声道

channel=str2double(get(handles.channel,'String'));

%将指定声道的信息存放到handles.y中

handles.y=temp(:,channel);

3.2.3 信号发生器

MATLAB有产生标准信号的函数，如sawtooth能够产生三角波或钜齿波，首先利用get函数获得波形soundtype，频率frequency，幅值amp和相位phase，然后是以下代码。

switch soundtype

    case 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 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');

end

if get(handles.add,'Value')==0.0

    handles.y=y;       %若没有勾选上“混迭”，则将生成的波形赋给handles.y

else                            %否则将生成的波形与原有波形叠加

    handles.y=handles.y+y;

end

3.3分析模块

由于MATLAB的绘图功能很强大，所以图形显示模块不用单独开发，可直接调用plot、axis等函数实现图形显示功能，故图形显示也将在分析模块中给出。

3.3.1 时域分析

2.1.2节给出时域分析中的过零检测算法流程，故这里不给出过零检测的代码。MATLAB提供了mean，std函数，能够方便地计算均值、标准差。下面是过零检测之后的代码，其中T为过零检测得到的周期(向量)，amp为过零检测得到的幅值(向量)，n为过零点数。

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);   %from，to即是界面中的“从第from点到第to点”

%计算并输出均值

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);

3.3.2 频域分析

频域分析需要作Fourier变换，MATLAB提供了fft函数，能够方便地实现快速Fourier变换算法。以下代码省去了从界面中获得from、to、Fs的部分，也省去了绘图后设置横、纵坐标轴的名称的部分。

%首先提取出待分析的样本，将其存入sample中

sample=handles.y(from:to);

%生成离散化的频率点，以采样频率作为离散化的间隔

f=linspace(0,Fs/2,(to-from+1)/2);

%对样本作快速Fourier变换，变换结果存入Y中

Y=fft(sample,to-from+1);

[C,I]=max(abs(Y));    %获得幅值最大的点及其所对应的下标值I

%则f(I)为最大的幅值所对应的频率，即信号频率的估计值

set(handles.foutt,'String',1/f(I)); %计算并输出周期的估计值

set(handles.foutfreq,'String',f(I));      %输出频率的估计值

Y=Y(1:(to-from+1)/2);                    %为与f对应，只取Y的前半部分

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);       %绘制功率谱曲线

4．运行实例与误差分析

为了分析软件的性能并比较时域分析与频域分析各自的优势，本章给出了两种分析方法的频率估计的比较，分析软件的在时域和频域的计算精度问题。

4.1频率估计

4.1.1 标准正弦信号的频率估计

用信号发生器生成标准正弦信号，然后分别进行时域分析与频域分析，得到的结果如图 4所示。从图中可以看出，时域分析的结果为f = 400.3702Hz，频域分析的结果为f = 417.959Hz，而标准信号的频率为400Hz，从而对于标准信号时域分析的精度远高于频域分析的精度。

 

图 4 标准正弦信号的频率估计

4.1.2 带噪声的正弦信号的频率估计

先成生幅值100的标准正弦信号，再将幅值50的白噪声信号与其混迭，对最终得到的信号进行时域分析与频域分析，结果如图 5所示，可以看出，时域分析的结果为f = 158.9498Hz，频域分析的结果为f = 200.391Hz，而标准信号的频率为200Hz，从而对于带噪声的正弦信号频域分析的精度远高于时域分析的精度。

图 5 带噪声的正弦信号的频率估计

4.2 结果分析与结论

在时域，频率估计是使用过零检测的方式计算出，从而对于带噪声的信号既容易造成“误判”，也容易造成“漏判”，且噪声信号越明显，“误判”与“漏判”的可能性越大。但在没有噪声或噪声很小时，时域分析对每个周期长度的检测是没有累积误差的，故随着样本容量的增大，估计的精度大大提高。

在频域，频率估计是通过找出幅值谱峰值点对应的频率求出。故不会有时域分析的问题。但频率离散化的误差及栅栏效应却是不可避免地带来误差，仅频率离散化的误差就大于Fs/2。

由实验结果及以上的分析可以得出结论：在作频率估计时，如果信号的噪声很小，采用时域分析的方法较好；如果信号的噪声较大，采用频域分析的方法较好。

5．总结与展望

5.1总结

本文给出了基于MATLAB的声音信号频谱分析仪的设计原理与实现方法，在原理部分，从时域和频域两个方面提供了信号分析所需要的算法流程及计算公式，在原理的最后还结合软件工程理论给出了软件的模块划分，这样在基于此设计原理的基础上可以用任何平台任何语言进行软件开发。在实现方法上，结合软件的界面和具体的代码讲述了整个软件编码实现的原理。最后结合一个运行实例比较了时域分析与频域分析计算频率的异同之处，并分析了误差的原因。

尽管MATLAB有强大的数学函数库，使得编程时间大大缩短，但MATLAB有它固有的缺陷，如运行速度太慢，因为它是解释型语言，而且运行依赖了MATLAB软件，无法发布为商用软件，另外在控制用户输入上也比较难以实现。这些缺陷也导致了用MATLAB所开发的软件有这些缺陷。

5.2展望

对于该软件，还有许多可以改进的地方。MATLAB做界面时难以控以用户的错误输入，但可以通过添加适当的判断语句提示用户的错误输入。对于软件运行依赖于MATLAB的问题，可以通过MATLAB与VC混合编程的方式解决，VC是Windows平台下主要的应用程序开发环境之一，它能方便实现软件开发，开发的系统具有界面友好、执行速度快、易维护和升级等优点，故两者的结合能提高开发效率而不损软件性能。

参考文献

[1]          黄昱. 软件频谱仪设计. 华中科技大学硕士学位论文. 2004.5

[2]          肖旸，胡耀祖. 基于虚拟技术的频谱分析仪器方案. 微机发展，2003.6

[3]          张桂林，张烈平. 基于声卡和Matlab的虚拟信号发生器. 现代电子技术

[4]          张海藩. 软件工程. 北京：人民邮电出版社, 2002.3

 

源代码下载http://d.download.csdn.net/down/364891/sbtdkj1017