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

哈尔滨工程大学

实验报告

实验名称：离散时间滤波器设计

班级：电子信息工程4班

学号：

姓名：

实验时间：2016年10月31日18：30

成绩：________________________________

指导教师：栾晓明

实验室名称：数字信号处理实验室哈尔滨工程大学实验室与资产管理处制

实验七音频频谱分析仪设计与实现

一、实验原理

MATLAB 是一个数据分析和处理功能十分强大的工程实用软件，其数据采集工具箱为实现数据的输入和输出提供了十分方便的函数命令。本实验要求基于声卡和MTLAB 实现音频信号频谱分析仪的设计原理与实现，功能包括：

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

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

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

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

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

3、相位检测

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

频率、幅值和相位估计的流程如图1所示。 4、数字信号统计量估计 (1) 峰值P 的估计

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

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

(2)均值估计

i N

i y N

y E ∑==

1

)(

式中，N 为样本容量，下同。 (3) 均方值估计

()20

2

1

∑==

N

i i

y

N

y E

(4) 方差估计

∑=-=

N

i i

Y E y N

y D 0

2

))

((1

)(

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

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

5、频谱分析原理

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

(1)DFT 与FFT

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

()t t f y f y Y N i i i i ???

?????+?=∑∑-101

-N 0)2sin(j t)2cos(f ππ

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

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

1,,2,1,0,)(/21

-=?=?--∑N k t e y f k Y N ki j N i π

式中，Δ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

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','Chooseerrer');

end

if get(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([0 N -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('No wave exist! Please choose a input type!');

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

if from<1 | to-from<5;

msgbox('Error range!');

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

ifymax

ymax=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('No wave exist! Please choose a input type!');

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，从而对于带噪声的正弦信号频域分析的精度远高于时域分析的精度。

结果分析：

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

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

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

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

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

(3)声卡输入

由麦克风录下自己读课本的一段声音，进行频域分析，得到周期为无穷大，如图6

图6声卡输入频域分析

(4)WAV文件输入

打开同一文件夹下的名为”冲击声”的WAV文件，进行时域和频域分析，结果如图7

图7WAV输入时域及频域分析

(5)波形发生器产生其他波形

输入方式为波形发生器，验证产生其他波形，并进行时域及频域分析，图8，9，10分别为方波，三角波，锯齿波的时域及频域分析结果。

图8方波时域及频域分析

图9三角波时域及频域分析

图10锯齿波时域及频域分析

四、实验总结

有程序运行结果及对设计要求功能的一一验证可知该GUI设计基本满足要求，能对不同输入信号进行时域及频域相关的各种分析，功能及其强大，对于信号处理领域有重要的意义，同时对于时域和频域分析的比较可知，在作频率估计时，如果信号的噪声很小，采用时域分析的方法较好；如果信号的噪声较大，采用频域分析的方法较好。

本次实验难度较大，前期准备花了很多实验才得到了这样一个满意的结果，虽然过程及其艰辛，但是收获却是极其丰富，对于matlab这一强大的数学分析软件有了更加深入的了解，掌握了GUI设计的过程及方法，对于以后的学习有着极大的帮助。同时此次设计的频谱分析仪对于今后信号处理相关知识的理解也可起到一定的辅助作用。