基于连续周期信号时域采样数据的频谱分析程序设计.doc
《基于连续周期信号时域采样数据的频谱分析程序设计.doc》由会员分享，可在线阅读，更多相关《基于连续周期信号时域采样数据的频谱分析程序设计.doc（11页珍藏版）》请在冰点文库上搜索。
陕西理工学院毕业设计论文
题目基于连续周期信号时域采样数据的频谱分析程序设计
学生姓名黄冬冬
学号1110064069
所在学院物理与电信工程学院
专业班级电信1103班
指导教师龙姝明
完成地点博远楼A1109，C1207
2015年6月5日
基于连续周期信号时域采样数据的频谱分析程序设计
黄冬冬
（陕西理工学院物理与电信工程学院电子信息科学与技术专业,2011级03班陕西汉中723000）
指导教师：
龙姝明
[摘要]大量的物理系统中，时域信号都是连续带限信号，即含有有限多个频率成分的连续信号，在观测之前，变化规律都是未知。
要探索这些带限信号的时域演化规律，必须对信号观察测量大量数据，即对信号采样产生大量采样数据，通过对信号的正确采样和对采样数据的科学频谱分析获得采样数据的频谱并证明它就是原连续周期信号的频谱。
本文分析了用数值计算的方法实现连续周期时间信号的时域频谱分析，即采用快速傅里叶变换（FFT）方法来获取连续周期信号的频谱，并编写实用的连续带限信号频谱分析Mathematica程序，以实例展示程序的高效率和可靠性。
[关键词]连续带限信号；采样数据；频谱分析；FFT；Mathematica程序
TheProgrammingofFrequencySpectrumAnalysisonContinuousPeriodicSignalSamplingData
HuangDongdong
（Grade11,Class03,MajorElectronicInformationScienceandTechnology,SchoolofPhysicsandTelecommunicationEngineering，ShaanxiUniversityofTechnology,Hanzhong723000,Shaanxi）
Tutor:
LongShuming
AbstractAlotofphysicalsystem,thetimedomainsignaliscontinuousband-limitedsignal,namelythecontinuoussignalcontainingmultiplefrequencycomponentsco.,LTD.,beforetheobservation,changelawisunknown.Toexplorethetime-domainevolutionoftheband-limitedsignal,mustwatchthesignalmeasurementoflargeamountsofdata,namelythesignalsamplingproducesalargenumberofsamplingdata,bymeansofthecorrectsignalsamplingandthesamplingdataofspectrumanalysistoobtainthefrequencyspectrumofthesampleddatascientificallyandprovethatitistheoriginalcontinuousspectrumofperiodicsignal.Analyzedinthispaperthemethodofnumericalcalculationisusedtoimplementcontinuouscycletimesignalspectrumanalysisoftimedomain,whichusesthefastFouriertransform（FFT）methodtogetthecontinuousspectrumofperiodicsignal,andwriteusefulforband-limitedsignalspectrumanalysisMathematicaprogram,withexamplestoshowtheefficiencyandreliabilityoftheprogram.
KeywordsContinuousband-limitedsignal,Samplingdata,Spectrumanalysis,FFT,Mathematicaprogram
目录
1已知变化规律的周期信号频谱分析方法 1
1.1周期信号的分解 1
1.2傅里叶级数展开条件 2
1.3指数形式傅里叶级数 2
1.4周期信号傅里叶级数展开的技巧性方法 2
1.4.1平移方法 2
1.4.2求导方法 2
1.5已知变化规律的周期连续信号频谱识别方法 3
2未知变化规律的周期信号的频谱分析方法 3
2.1用FFT计算未知变化规律的周期信号的频谱的思路 3
2.2未知变化规律的周期连续信号频谱识别方法 4
3周期信号频谱分析的Mathematica程序设计思路 4
4周期信号频谱分析的Mathematica程序应用实例 4
结语 7
附录 8
在现实生活中对于信号进行频谱分析具有重要的意义。
通过对信号频谱的分析，可以得到信号的频率结构，了解信号的频率成分或系统的特征。
在此基础之上，可实现对信号的跟踪控制，从而实现对系统状态的早期预测，发现潜在的危险并诊断可能发生故障的原因，对系统参数进行识别及校正。
因此，频谱分析是揭示信号特征的重要方法，也是处理信号的重要手段[1]。
对于已知变化规律的周期信号做频谱分析，可采用傅里叶级数展开方法和时域采样频谱分析方法。
使用时域采样频谱分析方法时，对采样得到的数据进行快速傅里叶变换，然后利用快速傅里叶变换的数据来计算出信号的频谱的每个频率成分的振幅、频率、初相位。
对于未知变化规律的周期信号做频谱分析只有唯一一种方法，即利用信号的时域采样数据取FFT来得到信号的频谱的方法[2]。
在工程领域中，Mathematica软件是一种倍受程序开发人员青睐的语言，对于一些需要做大量数据运算处理的复杂应用Mathematica软件显得游刃有余[3]。
1已知变化规律的周期信号频谱分析方法
1.1周期信号的分解
设有周期信号，他的周期是T，角频率,它可分解为
（1-1）
式（1-1）中的系数，称为傅里叶系数[4]。
它可由式
（1-2）
求得。
为简便，式的积分区间取为或。
由式（1-1）可得傅里叶系数
n=0，1，2，…（1-3）
n=0，1，2，…（1-4）
式中T为函数的周期，为角频率。
由式（1-3）和式（1-4）可见，傅里叶系数和都是n的函数，其中是n的偶函数，即；而是n的奇函数，即有。
将式（1-1）中同频率项合并，可写成如下形式
（1-5）
式中
n=1，2，…（1-6）
如将式（1-6）的形式转化为式（1-7）的形式，他们系数之间的关系为
n=1，2，…（1-7）
由式（1-7）可见，是n的偶函数，即有；而是n的奇函数，即有。
傅里叶系数的这些重要性质是很有用的。
1.2傅里叶级数展开条件
周期信号应满足狄里赫利条件，即：
（1）在其一个周期内绝对可积；
（2）在其一个周期内只有有限个有限的不连续点；
（3）在其一个周期内只有有限个极大值和极小值。
注意：
条件
（1）为充分条件但不是必要条件；
条件
（2）（3）是必要条件但不是充分条件[5]。
1.3指数形式傅里叶级数
连续时间周期信号可以用指数形式傅里叶级数表示为
（1-8）
n=±1两项的基波频率为f0，两项合起来称为信号的基波分量；
n=±2的基波频率为2f0，两项合起来称为信号的2次谐波分量；
n=±N的基波频率为Nf0，两项合起来称为信号的N次谐波分量。
物理含义：
周期信号可以分解成不同频率虚指数之和[6]。
1.4周期信号傅里叶级数展开的技巧性方法
1.4.1平移方法
有些周期信号的周期规律不明显或周期变化计算过于复杂，则可通过左右平移和增加（减小）直流增量来改变信号位置，由此可使信号由变化规律不明显信号变成奇信号或者偶信号，然后进行傅里叶级数展开时就会比较快。
设已知变化规律的周期信号为，那么对该周期信号进行左右平移即，或是上下平移即进行变换的话，则可以将一个变化规律不明显或者计算过程过于复杂的信号变成一个奇信号或偶信号。
变为奇信号以后，变为偶信号后，那么对于计算傅里叶级数的系数，就会简单许多。
1.4.2求导方法
有些周期信号的周期规律不明显或周期变化计算过于复杂，则可对其进行求导（求导只是对其一个周期求导，因为信号的所有信息在一个周期均有完全的反应），这样就是将复杂的信号分析改变成一个周期内的函数进行分析，可以减少计算过程，若经过一次求导以后函数依旧复杂，可再次对其进行求导。
给出一个周期的表达式，，周期信号表达式为
（1-9）
对其求一次导则有
（1-10）
转换可得：
（1-11）
对其求二次导则有
（1-12）
转换可得：
（1-13）
1.5已知变化规律的周期连续信号频谱识别方法
通过上面分析我们已经知道，对于这类信号要获得其频谱，我们只需将此周期连续信号做傅里叶级数展开[7]。
设f（t）为一个周期连续信号，T为其周期，f=1/T是信号的基波频率，为信号的基
波角频率，则f（t）可以表示为，（n=0,1,2………）。
对f（t）进行傅里叶级数展开。
（1-14）
所谓周期信号的频谱，就是指周期信号可以表达成一系列正弦分量的叠加，每一个正弦成分称为其周期信号的一个频率成分，其中需要用三个实数来描述，即振幅、频率、相位。
只要得到了每一个频率成分的振幅、频率、相位，那么叠加起来就是该周期函数的函数表达式，同时也可以得到任意时区上的波形[8]。
（1-15）
通过上面推导即可计算已知变化规律的周期连续信号的频谱和振幅，从而画出相应信号f（t）的幅频特性杆状图。
同时也可以对已知变化规律的周期信号分析其频谱时，也可以采用时域采用方法，即对周期信号的一个周期的时区上进行采样，利用采样数据的快速傅里叶变换再来计算每个频率成分的振幅、频率、初相位，也可以得到周期信号的频谱。
2未知变化规律的周期信号的频谱分析方法
对于未知变化规律的周期信号，首先可以确定信号是周期信号，虽然其周期和变化规律均不知道，但是其周期是客观存在的，可以使用设备对其进行测量的。
因为信号的周期不知道，所以不能采用傅里叶级数展开的方法，只能利用信号时域采样的数据取FFT来得信号其频谱的方法。
由于周期信号理论上用无穷个频率成分，进行时域采样就不能满足采样定理，即周期信号的采样频率必须大于信号最高频率的二倍，表达式为。
但是周期信号的高频成分的振幅都非常小，从级数层面上讲是可以忽略不计的，这样就可以把周期信号的最高频率堪称是有限的，只有这样才能使用时域采样的方法进行频谱分析[9]。
2.1用FFT计算未知变化规律的周期信号的频谱的思路
利用FFT对连续时间未知变化规律的周期信号进行频谱分析[10]其实是一个对信号进行逐级近似的过程：
（1）选取尽可能小的采样周期与尽可能宽的采样时区；
（2）选取尽可能多的采样数据点，对信号进行采样。
（3）画采样数据的波形图，初步判断它的周期，并用数学手段找出其精确周期。
（4）按照第三步得到的信号周期计算出一个周期的采样数据量，并截出这些数据。
（5）对选取的数据做快速傅里叶变换。
利用采样数据的快速傅里叶变换的结果计算信号的频率成分的振幅、频率、初相位，从而得到周期信号的频谱。
2.2未知变化规律的周期连续信号频谱识别方法
未知变化规律的周期连续信号通常是通过对信号采样，然后分析采样信号的频谱来获得连续信号的频谱，这种做法首先要满足采样定理。
为了使实际信号在采样后能够不失真的再现，采样频率必须大于信号最高频率的两倍。
其次，为获得连续信号的频谱，要对采样数据进行离散傅里叶变换（DFT）来获得离散信号的频谱[11]。
我们假设周期连续信号采样得到的序列为x（n），设x（n）为有限长序列，长度为N，但它取自周期信号，我们可以认为在采样参数选择恰当时，对它进行周期延拓所得周期序列的最小周期与被采样信号的周期是一致的。
通过采样定理的学习我们得知采样信号的频谱是被采样信号频谱的周期延拓，所以我们只需要对一个采样周期进行研究便可得到信号的全部信息。
现在我们取定采样序列x（n），长度为N（称为一个主值区间）对它进行频谱分析即可。
对x（n）取DFT，计算过程如下：
具体计算步骤如下：
（2-1）
（2-2）
式中是待采样得周期连续信号，T是采样周期，是采样后的有限长离散序列，与是等同的。
由上式可知有限长离散信号频谱的采样是由其连续频谱经过采样得到的。
由此得DTFT是计算离散序列的频谱函数的，而DFT则是计算离散序列频谱的采样的。
3周期信号频谱分析的Mathematica程序设计思路
首先定义一个周期信号，给出采样时间和采样点，画出信号在一个周期的波形图（周期信号的所有信息在一个周期内均有表达），然后对这一个周期进行频谱分析[11]。
（1）选择合适的采样间隔和合理的采样数据总量，对信号进行采样；
（2）对信号进行采样，记录测量的时间点和测量的值；
（3）计算信号的周期，采样数据话波形图，观察获取周期；
（4）截取一个周期的采样数据量，做FFT；
（5）利用变换的结果，计算前m个频率成分的振幅、频率、相位。
并筛选大振幅的频率成；
（6）验证：
利用保留的频率成分重新在时域里构造一个多频率正弦信号叠加的函数关系，然后画出波形图，再与采样数据画出的波形图进行比对。
如果符合，那么频谱分析结论就是可靠的。
4周期信号频谱分析的Mathematica程序应用实例
利用Mathematica软件对周期三角波信号进行频谱分析。
利用Mathematica软件编程，程序见附录。
运行程序可得出以下结论：
一个周期内的三角波形图如图1：
图1一个周期内的三角波形图
频率、振幅、相位分别如表1：
表1周期信号的前三个频率成分的频率、振幅、相位
名称
fo
f1
f2
频率
0
1000
3000
振幅
0.5
0.4054198
0.0451652
相位
0
3.14159
3.14159
重构信号y（t）=0.5+0.405418Cos[3.14159+2000πt]+0.0451652Cos[3.14159+6000πt]
三角波和正弦波在同一坐标系里的比较图如图2：
图2三角波和正弦波在同一坐标系里的比较图
最大误差绝对值为：
error=0.0494167
结语
以上对周期三角波信号的频谱进行了分析，在分析过程中需要对信号进行采样和数据计算，因此会出现一定误差，但只要合理选择分析参数可以使误差保持在工程允许的范围之内。
对于随机信号，由于是无限大能量的功率信号，它不满足傅里叶变换条件，而且也不存在解析表达式，因此就不能够应用确知信号的频谱计算方法去分析随机信号的频谱，但可以应用时域采样频谱分析方法对此信号进行频谱分析，通过上机编程可以解决该问题。
致谢
龙姝明老师在毕业设计中的指导和帮助，以及我多次找他他都不厌其烦地给解答我的各种困惑，细心、耐心地给我提出论文和程序里的多个问题，让我下来再继续修改和完善，使得我能够顺利地完成此次毕业设计。
通过本次Mathematica频谱分析毕业设计，让我了解了关于Mathematica软件在频谱分析方面的应用，又一次学习了Mathematica软件的使用和程序的设计，Mathematica的分析计算过程使我更加深入地了解了时域频谱分析的过程，使得我对频谱分析的理解更加深了一步。
Mathematica拥有强大的图形功能和高精度的数值计算功能，在生产和研究中起着重要的作用。
这次毕业设计，使我懂得了理论与实际相结合是很重要的，只有理论知识是远远不够的，只有把所学的理论知识与实践相结合起来，从理论中得出结论，才能真正为社会服务，从而提高自己的实际动手能力和独立思考的能力。
这次课程设计终于顺利完成，在设计中遇到的运行和调试问题，最后在老师的耐心指导下，终于游逆而解。
在以后的学习过程中我要不段的学习，不断丰富自己的知识。
参考文献
[1]倪铭.浅谈频谱分析[J].淮南职业技术学院学报.2007,01:
50-52.
[2]WeiYang,ZhiqinZhao,ZaipingNie.FastFourierTransformMultilevelFastMultipoleAlgorithminRoughOceanSurfaceScattering[J].Electromagnetics,2009,297-309.
[3]尹继武,龙姝明.基于Mathematica程序的正弦阶梯波频谱分析[J].四川师范大学学报（自然科学版）,2006,03:
376-378.
[4]吴大正,杨林耀,张永瑞,等.信号与线性系统分析（第4版）[M].北京:
高等教育出版社,2011.120-121.
[5]张睿.“信号与系统”辅导课的教学研究与实践[J].合肥工业大学学报（社会科学版）,2010,03:
156-158.
[6]张锐,袁丽英.信号与系统课程中对周期信号频谱分析的理解[J].高等函授学报（自然科学版）,2012,02:
30-32.
[7]刘志松.基于小波分析的信号去噪方法[J].浙江海洋学院学报（自然科学版）,2011,02:
150-154.
[8]杨宇,贾永兴.“信号与系统”周期信号频谱的教学分析[J].中国电力教育,2013,32:
146-149.
[9]龙姝明,朱杰武,孙彦清,等.数学物理方法&Mathematica[M].西安:
陕西人民教育出版社,2002.61-64.
[10]ApplicationofBayesianapproachtohydrologicalfrequencyanalysis[J].ScienceChina（TechnologicalSciences）,2011,05:
1183-1192.
[11]李昌利,沈玉利.“信号与系统”课程教学中的几点思考[J].高教论坛,2008,03:
122-123.
[12]孙云龙,张卫东.深入Mathematica编程[J].淄博学院学报（自然科学与工程版）,2000,02:
29-31.
附录
T=0.001;
Δ[t_]=（2t/T）（UnitStep[t]-UnitStep[t-T/2]）+（2（T-t）/T）（UnitStep[t-T/2]-UnitStep[t-T]）;
fs=10^5;Ts=1/fs;
L=117;
ts=Range[0,L-1]*Ts;
x=Δ[Mod[ts,T]];
L=Length[x];
ts=Range[0,L-1]*Ts;
data=Transpose[{ts,x}];
ps=ListPlot[data,Joined->True]
p=Flatten[Position[x,0.0]];
If[Length[p]<2,Print["88888"];Exit[]];
n=p[[2]]-p[[1]];
xn=Take[x,{1,n}];
X=（2/n）Fourier[xn,FourierParameters->{1,-1}];
X[[1]]/=2;
A=Abs[X[[1;;n/2]]]//Chop;
Am=Max[A];thr=0.05;
nh=Range[0,n/2-1];
F=（nh/n）fs;
φ=Arg[X[[1;;n/2]]];
sp0=Transpose[{A,F,φ}];
sp=Select[sp0,#[[1]]/Am>=thr&];
Print["信号频率成分：
{{A1,F1,φ1},...}=",sp];
Lf=Length[sp];
y[t_]=Sum[sp[[k,1]]Cos[2πsp[[k,2]]t+sp[[k,3]]],{k,1,Lf}];
Print["重构信号y（t）=",y[t]];
py=Plot[y[t],{t,0,（L-1）*Ts},PlotStyle->Red,PlotRange->{All,{-0.5,2}}];
Show[py,ps]
yn=y[ts];
error=Max[Abs[x-yn]]；
Print["最大误差error=",error];
第7页共8页
升级会员 