基于MATLAB的频谱分析仪设计.docx
- 文档编号:4750223
- 上传时间:2022-12-08
- 格式:DOCX
- 页数:25
- 大小:528.69KB
基于MATLAB的频谱分析仪设计.docx
《基于MATLAB的频谱分析仪设计.docx》由会员分享,可在线阅读,更多相关《基于MATLAB的频谱分析仪设计.docx(25页珍藏版)》请在冰豆网上搜索。
基于MATLAB的频谱分析仪设计
基于MATLAB的信号频谱分析仪的实现
一、概述
信号处理几乎涉及到所有的工程技术领域,而频谱分析又是信号处理中一个非常重要的分析手段。
一般的频谱分析都依靠传统频谱分析仪来完成,价格昂贵,体积庞大,不便于工程技术人员的携带。
虚拟频谱分析仪改变了原有频谱分析仪的整体设计思路,用软件代替了硬件,使工程技术人员可以用一部笔记本电脑到现场就可轻松完成信号的采集、处理及频谱分析。
在工程领域中,MATLAB是一种倍受程序开发人员青睐的语言,对于一些需要做大量数据运算处理的复杂应用以及某些复杂的频谱分析算法MATLAB显得游刃有余。
本文将重点介绍基于MATLAB的虚拟频谱分析仪的设计。
本文设计的虚拟频谱分析仪的功能包括:
(1) 音频信号信号输入。
输入的途径包括从声卡输入、从WAV文件输入、从信号发生器输入;
(2) 信号波形分析。
包括幅值、频率、周期、相位的估计,并计算统计量的峰值、均值、均方值和方差等信息;
(3) 信号频谱分析。
频率、周期的估计,图形显示幅值谱、相位谱和功率谱等信息的曲线。
二、实验原理
2.1时域抽样定理
时域抽样定理给出了连续信号抽样过程中信号不失真的约束条件:
对于基带信号,信号抽样频率大于等于2倍的信号最高频率,即。
时域抽样是把连续信号变成适于数字系统处理的离散信号。
对连续信号以间隔T抽样,则可得到的离散序列为。
图1连续信号抽样的离散序列
若,则信号与的频谱之间存在:
其中:
的频谱为,的频谱为。
可见,信号时域抽样导致信号频谱的周期化。
(rad/s)为抽样角频率,为抽样频率。
数字角频率Ω与模拟角频率ω的关系为:
Ω=ωT。
2.2离散傅立叶变换(DFT)
有限长序列
的离散傅立叶变换(DFT)为
逆变换为
2.3快速傅立叶变换(FFT)
在各种信号序列中,有限长序列占重要地位。
对有限长序列可以利用离散傅立叶变换(DFT)进行分析。
DFT不但可以很好的反映序列的频谱特性,而且易于用快速算法(FFT)在计算机上进行分析。
有限长序列的DFT是其z变换在单位圆上的等距离采样,或者说是序列傅立叶的等距离采样,因此可以用于序列的谱分析。
FFT是DFT的一种快速算法,它是对变换式进行一次次分解,使其成为若干小数据点的组合,从而减少运算量。
MATLAB为计算数据的离散快速傅立叶变换,提供了一系列丰富的数学函数,主要有Fft、Ifft、Fft2、Ifft2,Fftn、ifftn和Fftshift、Ifftshift等。
当所处理的数据的长度为2的幂次时,采用基-2算法进行计算,计算速度会显著增加。
所以,要尽可能使所要处理的数据长度为2的幂次或者用添零的方式来添补数据使之成为2的幂次。
Fft函数调用方式:
Y=fft(X);
Y=fft(X,N);
Y=fft(X,[],dim)或Y=fft(X,N,dim)。
函数Ifft的参数应用与函数Fft完全相同。
2.4波形分析原理
2.4.1 信号频率、幅值和相位估计
(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过零检测法估算信号的周期、幅值和相位流程图
2.4.2 数字信号统计量估计
(1) 峰值P的估计
在样本数据x中找出最大值与最小值,其差值为双峰值,双峰值的一半即为峰值。
P=[max(yi)-min(yi)]/2
(2)均值估计
式中,N为样本容量,下同。
(3) 均方值估计
(4) 方差估计
2.5频谱分析原理
时域分析只能反映信号的幅值随时间的变化情况,除单频率分量的简单波形外,很难明确提示信号的频率组成和各频率分量大小,而频谱分析能很好的解决此问题。
由于从频域能获得的主要是频率信息,所以本节主要介绍频率(周期)的估计与频谱图的生成。
2.5.1 频率、周期的估计
对于Y(kΔf),如果当kΔf = f时,Y(kΔf)取最大值,则f为频率的估计值,由于采样间隔的误差,f也存在误差,其误差最大为Δf /2。
周期T=1/f。
从原理上可以看出,如果在标准信号中混有噪声,用上述方法仍能够精确地估计出原标准信号的频率和周期,这个将在下一章做出验证
2.5.2 频谱图
为了直观地表示信号的频率特性,工程上常常将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]
三、程序设计及编写
采样频率Fs与采样点数N是声音信号输入时共同需要作用的参数,故将其独立出来。
下面为别介绍三种输入方式的实现。
3.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 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.3 信号发生器输入
MATLAB有产生标准信号的函数,如sawtooth能够产生三角波或钜齿波,首先利用get函数获得波形soundtype,频率frequency,幅值amp和相位phase,然后是以下代码。
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; %若没有勾选上“混迭”,则将生成的波形赋给handles.y
else %否则将生成的波形与原有波形叠加
handles.y=handles.y+y;
end
3.4 时域分析
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.5 频域分析
频域分析需要作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.1标准正弦信号的频率估计
用信号发生器生成标准正弦信号,然后分别进行时域分析与频域分析,得到的结果如图 4所示。
从图中可以看出,时域分析的结果为f =400Hz,频域分析的结果为f =400.2Hz,而标准信号的频率为400Hz,从而对于标准信号时域分析的精度远高于频域分析的精度。
图3标准正弦波分析
4.2 非标准正弦信号的频率估计
先成生幅值600的标准正弦信号,再将幅值300的白噪声信号与其混迭,对最终得到的信号进行时域分析与频域分析,结果如图 5所示,可以看出,时域分析的结果为f =388.68Hz,频域分析的结果为f =400.21Hz,而标准信号的频率为400Hz,从而对于带噪声的正弦信号频域分析的精度远高于时域分析的精度。
图3非标准正弦波分析
4.3小结
过零检测的方式对于带噪声的信号既容易造成“误判”,也容易造成“漏判”,且噪声信号越明显,“误判”与“漏判”的可能性越大。
但在没有噪声或噪声很小时,时域分析对每个周期长度的检测是没有累积误差的,故随着样本容量的增大,估计的精度大大提高。
但把信号进行傅里叶变换后,频率估计是通过找出幅值谱峰值点对应的频率求出。
不存在零点误判的问题。
但频率离散化的误差及栅栏效应却是不可避免地带来误差。
因此,在作频率估计时,如果没有干扰信号或干扰信号很小时,采用时域分析的方法较好;如果干扰信号较大,采用频域分析的方法较好。
附录资料:
matlab画二次曲面
一、螺旋线
1.静态螺旋线
a=0:
0.1:
20*pi;
h=plot3(a.*cos(a),a.*sin(a),2.*a,'b','linewidth',2);
axis([-50,50,-50,50,0,150]);
gridon
set(h,'erasemode','none','markersize',22);
xlabel('x轴');ylabel('y轴');zlabel('z轴');
title('静态螺旋线');
2.动态螺旋线
t=0:
0.1:
10*pi;
i=1;
h=plot3(sin(t(i)),cos(t(i)),t(i),'*','erasemode','none');
gridon
axis([-22-22035])
fori=2:
length(t)
set(h,'xdata',sin(t(i)),'ydata',cos(t(i)),'zdata',t(i));
drawnow
pause(0.01)
end
title('动态螺旋线');
(图略)
3.圆柱螺旋线
t=0:
0.1:
10*pi;
x=r.*cos(t);
y=r.*sin(t);
z=t;
plot3(x,y,z,'h','linewidth',2);
gridon
axis('square')
xlabel('x轴');ylabel('y轴');zlabel('z轴');
title('圆柱螺旋线')
二、旋转抛物面
b=0:
0.2:
2*pi;
[X,Y]=meshgrid(-6:
0.1:
6);
Z=(X.^2+Y.^2)./4;
meshc(X,Y,Z);
axis('square')
xlabel('x轴');ylabel('y轴');zlabel('z轴');
shadingflat;
title('旋转抛物面')
或直接用:
ezsurfc('(X.^2+Y.^2)./4')
三、椭圆柱面
loadclown
ezsurf('(2*cos(u))','4*sin(u)','v',[0,2*pi,0,2*pi])
view(-105,40)%视角处理
shadinginterp%灯光处理
colormap(map)%颜色处理
gridon%添加网格线
axisequal%使x,y轴比例一致
xlabel('x轴');ylabel('y轴');zlabel('z轴');
shadingflat;
title('椭圆柱面')%添加标题
四、椭圆抛物面
b=0:
0.2:
2*pi;
[X,Y]=meshgrid(-6:
0.1:
6);
Z=X.^2./9+Y.^2./4;
meshc(X,Y,Z);
axis('square')
xlabel('x轴');ylabel('y轴');zlabel('z轴');
shadingflat;
title('椭圆抛物面')
或直接用:
ezsurfc('X.^2./9+Y.^2./4')
五、'双叶双曲面
ezsurf('8*tan(u)*cos(v)','8.*tan(u)*sin(v)','2.*sec(u)',[-pi./2,3*pi./2,0,2*pi])
axisequal
gridon
axissquare
xlabel('x轴');ylabel('y轴');zlabel('z轴');
shadingflat;
title('双叶双曲面')
六、双曲柱面
loadclown
ezsurf('2*sec(u)','2*tan(u)','v',[-pi/2,pi/2,-3*pi,3*pi])
holdon%在原来的图上继续作图
ezsurf('2*sec(u)','2*tan(u)','v',[pi/2,3*pi/2,-3*pi,3*pi])
colormap(map)
shadinginterp
view(-15,30)
axisequal
gridon
axisequal
xlabel('x轴');ylabel('y轴');zlabel('z轴');
shadingflat;
title('双曲柱面')
七、双曲抛物面(马鞍面)
[X,Y]=meshgrid(-7:
0.1:
7);
Z=X.^2./8-Y.^2./6;
meshc(X,Y,Z);
view(85,20)
axis('square')
xlabel('x轴');ylabel('y轴');zlabel('z轴');
shadingflat;
title('双曲抛物面')
或直接用:
ezsurfc('X.^2./8-Y.^2./6')
八、抛物柱面
[X,Y]=meshgrid(-7:
0.1:
7);
Z=Y.^2./8;
h=mesh(Z);
rotate(h,[101],180)%旋转处理
%axis([-8,8,-8,8,-2,6]);
axis('square')
xlabel('x轴');ylabel('y轴');zlabel('z轴');
shadingflat;
title('抛物柱面')
或直接用:
ezsurfc('Y.^2./8')
九、环面
ezmesh('(5+2*cos(u))*cos(v)','(5+2*cos(u))*sin(v)','2*sin(u)',[0,2*pi,0,2*pi])
axisequal
gridon
xlabel('x轴');ylabel('y轴');zlabel('z轴');
shadingflat;
title('环面')
十、椭球
ezsurfc('(5*cos(u))*sin(v)','(3*sin(u))*sin(v)','4*cos(v)',[0,2*pi,0,2*pi])
axisequal
gridon
xlabel('x轴');ylabel('y轴');zlabel('z轴');
shadingflat;
title('椭球')
十一、单叶双曲面
ezsurf('4*sec(u)*cos(v)','2.*sec(u)*sin(v)','3.*tan(u)',[-pi./2,pi./2,0,2*pi])
axisequal
gridon
xlabel('x轴');ylabel('y轴');zlabel('z轴');
shadingflat;
title('单叶双曲面')
十二、旋转单叶双曲面
loadclown
ezsurf('8*sec(u)*cos(v)','8.*sec(u)*sin(v)','2.*tan(u)',[-pi./2,pi./2,0,2*pi])
colormap(map)
view(-175,30)
%alpha(.2)%透明处理
axisequal
gridon
axissquare
xlabel('x轴');ylabel('y轴');zlabel('z轴');
shadingflat;
title('旋转单叶双曲面')
十三、圆柱面
subplot(1,2,1)
ezsurf('(2*cos(u))','2*sin(u)','v',[0,2*pi,0,2*pi])
gridon
shadinginterp
axisequal
xlabel('x轴');ylabel('y轴');zlabel('z轴');
title('圆柱面')
subplot(1,2,2)
cylinder(30)
shadinginterp
axissquare
title('调用cylinder函数所得圆柱面')
十四、二次锥面
clc,clear;
P=[1,0,0;
0,cos(45*pi/180),sin(45*pi/180);
0,-sin(45*pi/180),cos(45*pi/180)];
fork2=1:
31
fork1=1:
31
x(k1,k2)=(k2-1)*cos((k1-1)*12*pi/180);
y(k1,k2)=(k2-1)*sin((k1-1)*12*pi/180);
z(k1,k2)=sqrt(x(k1,k2)^2+y(k1,k2)^2);
Vxyz=P*[x(k1,k2),y(k1,k2),z(k1,k2)]';
x1(k1,k2)=Vxyz
(1);
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 基于 MATLAB 频谱 分析 设计