9种谱校正方法
9种谱校正方法及matlab代码

9种谱校正方法及matlab 程序代码采样间隔归一化成1T ∆=,采样长度为N .这样FFT 离散谱线为0,1)i X i N =-(,相应的频率分辨率2/(1/)N f N ωπ∆=∆=. 设FFT 离散谱线局部极高谱线为m (为了数学上简洁,假定从0开始,注意在MATLAB 环境下数组实际操作的是从1开始),记频偏量δωδω=∆. 我们需要使用谱线m 和与之相邻一条次高谱线,记这连续两条谱线中左边一条序号为M (当次高谱线在m 左侧时1M m =-,反之M m =).下面列出若干算法的δ计算公式1. 加矩形窗的精确谱校正[1]i i i X U jV =+111()sin()()cos()M M M M opt M MV V M U U M K U U ωω+++-∆+-∆=- 1211cos()sin()cos()sin()opt M M opt M M K M Z V U M K M Z V U M ωωωωωω++-∆⎡⎤=+⎢⎥∆⎣⎦-∆+∆⎡⎤=+⎢⎥∆+∆⎣⎦ 2121cos()cos()()Z M Z M M m Z Z ωωωδ∆+∆-∆=+-- 2. 加矩形窗情形,采用解析单频模型的幅值比校正[1, 2]11||()||||M M M X M m X X δ++=+-+ 3. 加汉宁窗情形,采用解析单频模型的幅值比校正[1, 2]112||||()||||M M M M X X M m X X δ++-=+-+ 4. 加矩形窗情形,采用解析单频模型的复比值校正[3]11Re ()M M M X M m X X δ++⎛⎫=+- ⎪-⎝⎭5. 加汉宁窗情形,采用解析单频模型的复比值校正[3]112()M M M MX X M m X X δ+++=+-- 6. 加矩形窗情形,采用解析单频模型的复合复比值校正[3]11Re ()M m M M X M m X X δ++⎛⎫=+- ⎪-⎝⎭ 11m R m mX X X δ++=-,1111m m L m m m m X X X X X X δ---=-=-- 0.5)0.5)m L m R δδδδδ=-++((7. 加汉宁窗情形,采用解析单频模型的复合复比值校正[3]112Re ()M M m M M X X M m X X δ++⎛⎫+=+-⎪-⎝⎭ 112m m R m mX X X X δ+++=-,1111221m m m m L m m m m X X X X X X X X δ----++=-=-- 0.5)0.5)m L m R δδδδδ=-++((8. 加矩形窗,Quin 校正[4]11Re()Re(),Re()Re()m m L R m m X X X X αα-+== 11L R L R L Rααδδαα==---,, 00 R R L R δδδδδ>>⎧=⎨⎩当且其它9. 加汉宁窗,Quin 校正[4]11Re()Re(),Re()Re()m m L R m m X X X X αα-+== 212111L RL R L Rααδδαα++==---,,00 R R L R δδδδδ>>⎧=⎨⎩当且其它References1. Schoukens, J., R. Pintelon,H. Van Hamme. The interpolated fast Fourier transform: Acomparative study . IEEE Transactions on Instrumentation and Measurement. 1992, 41(2):226-232.2. 谢明,丁康. 频谱分析的校正方法.振动工程学报. 1994, 7(2): 172-179.3. 陈奎孚, 王建立,张森文. 频谱校正的复比值法.振动工程学报(已投). 2007.4. Quinn, B.G. Estimating frequency by interpolation using Fourier coefficients.IEEETransactions on Signal Processing. 1994, 42(5): 1264-1268.%========================这是调用调试==================DT=1;N=1024;PHI=pi/3;Ampl=1;CiR=11.9; %cycles in recordFreq=CiR/(DT*N); %frequencyTV=[0:N-1];DatVec=Ampl*cos(Freq*TV*2*pi+PHI);FV=fft(DatVec);figuresubplot(2,1,1);plot(TV,DatV ec);subplot(2,1,2);plot(abs(FV(1:round(N/2.56))));grid on[MV,MI]=max(abs(FV));%加矩形窗的解析校正--1FreqShift=SpecCorr(FV,MI,N,1);%加矩形窗的解析单频模型校正--2FreqShift=SpecCorr(FV,MI,N,2);%加汉宁窗的解析单频模型校正--3HanDat=DatVec.*hanning(N,'periodic')';FV=fft(HanDat);FreqShift=SpecCorr(FV,MI,N,3);%加矩形窗的解析单频模型校正+复比值法--4FV=fft(DatVec);FreqShift=SpecCorr(FV,MI,N,4);%加汉宁窗的解析单频模型校正+复比值法--5HanDat=DatVec.*hanning(N,'periodic')';FV=fft(HanDat);FreqShift=SpecCorr(FV,MI,N,5);%加矩形窗的解析单频模型校正+复比值法+左右平均--6FV=fft(DatVec);FreqShift=SpecCorr(FV,MI,N,6);%加汉宁窗的解析单频模型校正+复比值法+左右平均--7HanDat=DatVec.*hanning(N,'periodic')';FV=fft(HanDat);FreqShift=SpecCorr(FV,MI,N,7);%加矩形窗的解析单频模型校正+Quinn算法--8FV=fft(DatVec);FreqShift=SpecCorr(FV,MI,N,8);%加汉宁窗的解析单频模型校正+Quinn算法--9HanDat=DatVec.*hanning(N,'periodic')';FV=fft(HanDat);FreqShift=SpecCorr(FV,MI,N,9);===========这是子程序======================%spectrum correction assemble% the sampling interval is 1 s (or unitary)%Input: SpecVec--Discrte Fourier Spectrum from FFT% PeakIdx--the peak index, noting the matrix in MatLab start from 1% TL--the length (or the point number) of the FFT% method--correction method% output: PeakShift-- the corrected peak shifting from the peak in discrete% spectrumfunction PeakShift=SpecCorr(SpecVec,PeakIdx,TL,method)% picking up the second highest spectrum lineif(abs(SpecVec(PeakIdx-1))>abs(SpecVec(PeakIdx+1)))IP=[PeakIdx-1,PeakIdx];ShiftCorr=-1; %shift aligning with the PeakIdxelseIP=[PeakIdx,PeakIdx+1];ShiftCorr=0; %shift aligning with the PeakIdxendII=IP(1)-1; % noting that the index of a matrix in MATLAB starts from 1, not zero if(method==1) %an analyitic solution for rectangular windowU=real(SpecVec(IP));V=imag(SpecVec(IP));DW=2*pi/TL;KOPT=(sin(II*DW)*(V(2)-V(1))+cos(II*DW)*(U(2)-U(1)))/(U(2)-U(1));Z=V.*(KOPT-cos((IP-1)*DW))./(sin(DW*(IP-1)))+U;Tmp1=(Z(2)*cos(DW*(II+1))-Z(1)*cos(DW*II))/(Z(2)-Z(1));PeakPos=acos(Tmp1)/DW;PeakShift=PeakPos-(PeakIdx-1);elseif(method==2) %based on the analytical-single-tone model for rectangular window PeakShift=abs(SpecVec(IP(2)))/(abs(SpecVec(IP(2)))+abs(SpecVec(IP(1))));PeakShift=PeakShift+ShiftCorr; %shift aligning with the PeakIdxelseif(method==3) %based on the analytical-single-tone model for Hanning windowPeakShift=(2*abs(SpecVec(IP(2)))-abs(SpecV ec(IP(1))))/(abs(SpecV ec(IP(2)))+abs(SpecVec(IP(1 ))));PeakShift=PeakShift+ShiftCorr; %shift aligning with the PeakIdxelseif(method==4) %based on the analytical-single-tone model for rectangular window with complex correctionPeakShift=real(SpecVec(IP(2))/(SpecVec(IP(2))-SpecVec(IP(1))));PeakShift=PeakShift+ShiftCorr; %shift aligning with the PeakIdxelseif(method==5) %based on the analytical-single-tone model for Hanning window with complex correctionPeakShift=(2*SpecVec(IP(2))+SpecVec(IP(1)))/(SpecVec(IP(2))-SpecVec(IP(1)));PeakShift=real(PeakShift)+ShiftCorr; %shift aligning with the PeakIdx elseif(method==6) %based on the analytical-single-tone model for rectangular window with complex correction+averagePeakShift=real(SpecVec(IP(2))/(SpecVec(IP(2))-SpecVec(IP(1))));MaxPeakShift=PeakShift+ShiftCorr; %shift aligning with the PeakIdxLeftShift=real(SpecVec(PeakIdx)/(SpecVec(PeakIdx)-SpecVec(PeakIdx-1)))-1;RightShift=real(SpecVec(PeakIdx+1)/(SpecVec(PeakIdx+1)-SpecVec(PeakIdx)));%averagePeakShift=(0.5-MaxPeakShift)*LeftShift+(0.5+MaxPeakShift)*RightShift;elseif(method==7) %based on the analytical-single-tone model for Hanning window with complex correction+average,????PeakShift=(2*SpecVec(IP(2))+SpecVec(IP(1)))/(SpecVec(IP(2))-SpecVec(IP(1)));MaxPeakShift=real(PeakShift)+ShiftCorr; %shift aligning with the PeakIdxLeftShift=(2*SpecVec(PeakIdx)+SpecV ec(PeakIdx-1))/(SpecVec(PeakIdx)-SpecVec(PeakIdx-1))-1;RightShift=(2*SpecVec(PeakIdx+1)+SpecV ec(PeakIdx))/(SpecVec(PeakIdx+1)-SpecVec(PeakIdx ));%averagePeakShift=(0.5-MaxPeakShift)*LeftShift+(0.5+MaxPeakShift)*RightShift;elseif(method==8) %Quinn method for the rectangular windowa1 = real(SpecV ec(PeakIdx-1)/SpecVec(PeakIdx)); %lefta2 = real(SpecV ec(PeakIdx+1)/SpecVec(PeakIdx)); %rightLeftShift = a1/(1-a1);RightShift = -a2/(1-a2);if (LeftShift>0 & RightShift>0)PeakShift = RightShift;elsePeakShift = LeftShift;endelseif(method==9) %Quinn method for the Hanning windowa1 = real(SpecV ec(PeakIdx-1)/SpecVec(PeakIdx)); %lefta2 = real(SpecV ec(PeakIdx+1)/SpecVec(PeakIdx)); %rightLeftShift = (2*a1+1)/(1-a1);RightShift = -(2*a2+1)/(1-a2);if (LeftShift>0 & RightShift>0)PeakShift = RightShift;elsePeakShift = LeftShift;endendend。
第9章 紫外吸收光谱分析

讨论:
(1) 转动能级间的能量差Δ Ε r:0.005~0.050eV,跃迁 产生吸收光谱位于远红外区。远红外光谱或分子转动光谱; (2) 振动能级的能量差Δ Ε v约为:0.05~1eV,跃迁产 生的吸收光谱位于红外区,红外光谱或分子振动光谱; (3) 电子能级的能量差Δ Ε e较大1~20eV,电子跃迁产生 的吸收光谱在紫外—可见光区,紫外—可见光谱或分子的电 子光谱;
p → p*跃迁:红移; ;e
pp np
max(正己烷)
230 329
max(氯仿)
238 315
max(甲醇)
237 309
max(水)
243 305
溶剂的影响
苯
1
1:乙醚
酰
丙
2:水
酮
2
极性溶剂使精细结构 消失;
250 300
非极性 → 极性 n → p*跃迁:兰移; ;e p → p*跃迁:红移; ;e
(2)共轭烯烃中的 p → p*
p*
p*₃
p*
p p*
165nm 217nm p₂
(HOMO LVMO) p
p₁
p
max
共轭烯烃(不多于四个双键)p p*跃迁吸收峰位置可由伍德
沃德——菲泽 规则估算。 max= 基+nii
基-----是由非环或六环共轭二烯母体决定的基准值; 无环、非稠环二烯母体: max=217 nm
红移与蓝移
有机化合物的吸收谱带 常常因引入取代基或改变溶 剂使最大吸收波长λmax和吸 收强度发生变化:
λmax向长波方向移动称为 红移,向短波方向移动称为 蓝移 (或紫移)。吸收强度即 摩尔吸光系数ε增大或减小的 现象分别称为增色效应或减 色效应,如图所示。
化学杀虫剂中9种拟除虫菊酯检测方法

化学杀虫剂中9种拟除虫菊酯检测方法作者:孙长恩等来源:《江苏农业科学》2014年第04期摘要:采用HP-5MS毛细管柱和FID检测器,建立了同时测定农药杀虫剂中隐性成分四氟苯菊酯、烯丙菊酯、氯氟醚菊酯、联苯菊酯、高效氯氟氰菊酯、氯菊酯、高效氯氰菊酯、S-氰戊菊酯、溴氰菊酯等9种拟除虫菊酯的气相色谱分析方法。
9种拟除虫菊酯的标准偏差均小于0.03,变异系数小于3.0%,线性相关系数均大于0.999;9种拟除虫菊酯成分5个水平的添加回收率为95.5%~107.5%。
该法分离效果好、准确度高、重现性好且操作简单快速。
关键词:拟除虫菊酯;隐性成分;气相色谱中图分类号: S482.3+5 文献标志码: A 文章编号:1002-1302(2014)04-0273-02收稿日期:2013-08-13作者简介:孙长恩(1970—),男,江苏徐州人,硕士,高级工程师,从事农药产品质量检验及全分析工作。
Tel:(025)84470312;E-mail:**********************。
近年来,一些农药生产企业为了尽快使自制产品在同类产品中脱颖而出,获得更高效益,随意添加隐性成分现象比较严重。
生产商根据不同产品,针对不同使用目的,添加不同隐性成分,使其产品在防止同类害虫时,显示出比未添加隐性成分产品的“良好效果”[1]。
部分农药生产企业在农药中以增效剂名义加入另外一种农药的违规现象,特别是在全面禁止使用甲胺磷等高毒农药后出现一个较为突出的新问题,它可能成为现在或将来农产品安全重大隐患和农药管理的一个盲点[2]。
在农药中任意添加其他农药成分,不但违反了《农药管理条例》,更有可能对农业生产安全和病虫害治理产生重大隐患,甚至严重干扰我国农业重大生物灾害的防控行动。
为了维护广大消费者利益,保持农药杀虫剂可持续发展,近两年国家有关部门在实施农药国家监督抽查时对农药中隐性成分进行了检测,发现部分企业违规添加隐性成分。
光谱 校准 奇异值分解

光谱校准奇异值分解光谱校准是光谱分析中的一个关键环节,通过对光谱数据进行校准,可以提高光谱分析的准确性和可靠性。
其中,奇异值分解(Singular Value Decomposition, SVD)是常用的光谱校准方法之一。
本文将以简体中文介绍光谱校准和奇异值分解的原理、方法及其在光谱分析中的应用。
光谱是指在不同波长的电磁辐射下,物体所发射、吸收或散射的光的强度分布。
光谱分析是一种常用的分析手段,可以通过测量物体在不同波长下的光谱信息,来获取物质的结构、组成和性质等信息。
然而,光谱数据受到很多因素的影响,如仪器漂移、噪声、非线性等,这些影响会导致光谱数据的失真和偏差,从而影响光谱分析的准确性。
为了解决这些问题,光谱校准应运而生。
光谱校准是一种通过数学方法对光谱数据进行修正和优化的过程,主要目的是消除或减小仪器误差、噪声和其他影响因素对光谱数据的影响,从而提高光谱分析的准确性和可靠性。
奇异值分解是一种常用的矩阵分解方法,可以将一个矩阵分解为三个矩阵的乘积,其中一个矩阵是奇异值矩阵,另外两个矩阵是正交矩阵。
在光谱校准中,可以将光谱数据矩阵进行奇异值分解,通过对奇异值矩阵的处理,实现对光谱数据的校准和优化。
具体而言,光谱校准中的奇异值分解主要包括以下几个步骤:1.构建光谱数据矩阵:将采集到的光谱数据按照一定的格式组织成矩阵形式,其中每一行代表一个光谱样本,每一列代表一个波长点的光强值。
2.进行奇异值分解:对光谱数据矩阵进行奇异值分解,得到三个矩阵U、S和VT,其中U和VT是正交矩阵,S是奇异值矩阵。
3.选择合适的奇异值并修正:根据奇异值的大小来选择前几个重要的奇异值,并对其进行修正。
通常情况下,前几个奇异值代表了光谱数据中最重要的信息,因此可以选择这些奇异值进行修正。
4.重建光谱数据矩阵:通过修正后的奇异值和原始的正交矩阵U 和VT,重建光谱数据矩阵。
这样得到的重建矩阵可以更好地反映光谱数据的真实情况,消除了仪器漂移、噪声和非线性等因素的影响。
谱方法在计算流体力学中的应用研究

西北工业大学硕士学位论文谱方法在计算流体力学中的应用研究姓名:王建瑜申请学位级别:硕士专业:计算数学指导教师:欧阳洁200703013.3谱近似(a)右端函数(b)精确解图3.1右端函数与精确解3.3.1ChebyshevTau方法应用谱方法数值求解问题(3.1)式,首先要确定展开基函数。
根据第二章中关于基函数的介绍,选择Chebyshev多项式(同样可以选择Legendre多项式)作为基函数,其表达式为丸。
(x,Y)=乙(x)瓦∽其中乙(x)为m次Chebyshev多项式。
根据2.2.2小节中关于Chebyshev多项式的介绍,由于基函数不满足问题(3一1)式的边界条件,应用Tau方法来逼近。
为使基函数和权函数满足单位正交关系,选权函数为‰(x,力=最(x)6(y)其中删=石2再1獬铲∽冀此时有如下两种单位正交关系成立虬舭朋心y蚴={:!富剧硼£~鹕∽出={:!等“。
ⅣO,_y)=∑∑口射如(x,y)(3·8)将近似解(3-8)式代入问题(3—1)式的微分方程,并对余量进行加权积分眦等+等Mw)蛐=“m蒯_伽少,k,l=0,1,--.,N由权函数和基函数的正交性,得a£’o’+4眢’2’=厶,||},,=0,1,…,Ⅳ(3.9)(3-9)式中厶=£l£lj;f,盯(训),(x,y)dxdy,k,l=0,1,…,N而且(3.6)式和(3-7)式依然成立,应用这两个关系式代替函数关于两个空问变量的二阶偏导数的展开系数,则方程组(3-9)式等价于下列方程组去耋p2_k2h+去”p2_12ZP(Pk21y,PCP12‰=1‘1鼽似y№y)dxdy÷)%+■)%=rI∥甜(w),(而q;i:妊戢“pp+f为=l+2儡敢k,,=0,1,…,N(3一12)求解代数方程组(3-12)式,即可得(3.8)式所表示的近似解。
3.3.3结果分析首先研究右端函数基于经典Chebyshev多项式的展开。
九种常见的家谱格式

九种常见的家谱格式
1. 血簪家谱:
血簪家谱是最常见的家谱格式,通常按照祖孙正序排列,以男性为主,女性则标注婚配的夫家姓氏。
2. 家族谱:
家族谱是按照同一个祖宗的后代进行排列,包括所有后代的姓名、出
生日期、配偶、孩子等信息。
3. 地区谱:
地区谱是按照某一个地区的所有家族进行排列,其中包括这些家族的
发展历史、族谱、文化特点等内容。
4. 世系谱:
世系谱是按照一个家族的世系关系排列,主要是为了记录家族树的延
伸和发展情况。
5. 改元谱:
改元谱是按照某一位祖先的改元年号进行排列,以记录家族世代的延
续和变化。
6. 分支谱:
分支谱是指将同一个家族的后代按照分支进行排列,主要是为了记录家族在不同地区的分支发展情况。
7. 家谱图:
家谱图是将家族树的图像化表现,以方便人们更加清晰地了解家族发展的历程和世系关系。
8. 纪念录:
纪念录是一种较为简略的家谱,包括主要的家族成员信息以及重要事件、文化传承、家族口传历史等内容。
9. 家训:
家训是对家族的传统、文化和价值观进行总结和归纳,以指导后代继承和发扬家族优良传统和文化精神。
火花源原子发射光谱法测定线材中9种元素

A理但is骚-但字分册PTCA(PART B:CHEM.ANAL.)I工作简报DOI:10.11973lhjy-h\202l03012火花源原子发射光谱法测定线材中9种元素杨再军,刁正斌,王娟,羊绍松(攀枝花钢饥有限公司制造部.攀枝花617000)摘要:建立了火花源原子发射光谱法测定直径为8〜12mm的小规格线材中的碳、硅、猛、磷、硫、鎳l、珞、铜和筑等9种元素含量的方法。
试样经切割和打磨后.按规格及钢种选择对应的自制夹具固定.垂直于磨样机打磨试样的端面。
用随V型板专用的定位夹具将样晶定位在火花源原子发射光谱仪激发孔中心位置.使试样支架中心和电极中心对齐。
盖好顶盖.用氫气吹扫火花室5s.采用能量为0.2J、频率为100Hz激发光激发样品。
在优化的仪器工作条件下,试样中9种元素可在10min内完成测定。
采用块状光谱国家标准样品制作校准曲线.用校正公式消除了共存元素的干扰及由标准试样和实际试样形状差异带来的系统偏差。
结果显示:9种元素的质量分数均在一定范围内与其对应的响应值呈线性关系.检出限(3s)为1.5〜18.1“g•g1;用此方法分析实际样品.所得测定值与参考文献使用的其他方法的基本一致.测定值的相对标准偏差("=11)均在5%以内。
关键词:自制夹具;火花源原子发射光谱法;元素;小规格线材中图分类号:0657.31文献标志码:A文章编号:1001-4020(2021)03-0258-06热轧带肋钢筋等线材产品的国家标准要求分析其中的碳、硅、猛、磷、硫、镰、锯、铜和帆等9种元素'.因此.建立一种可以快速、准确地测定钢筋钢中化学成分的方法.对生产工艺及产品质量的控制具有重要意义。
通常采用国家标准方法湿法测定线材中的元素反。
对于直径为8〜12mm的线材,需要先用专用工具钻取装置钻取屑样.然后针对其中的碳、硫两元素采用高频感应燃烧-红外吸收光谱法进行测定.其余元素采用电感耦合等离子体原子发射光谱法(ICP-AES)进行测定'.该方法试样前处理流程复杂•分析周期长,难以满足线材产品大批量检验的需求。
等离子体质谱法测定Nb、Ta、Zr、Hf等9种元素量

等离子体质谱法测定Nb、Ta、Zr、Hf等9种元素量1 范围本方法规定了各类样品中Be、Cs、Hf、Li、Nb、Rb、Sr、Ta、Zr元素含量的测定方法。
本方法适用于矿物、水系沉积物、土壤和岩石等类型样品中以上各元素量的测定。
本方法检出限(3S)及测定范围见表1。
表1 方法检出限及测定范围*元素方法检出限(3S)μg·g-1测定范围μg·g-1元素方法检出限(3S)μg·g-1测定范围μg·g-1Be 0.006 0.02~200 Nb 0.03 0.1~5000Cs 0.003 0.01~200 Rb 0.1 0.3~5000 Hf 0.015 0.05~500 Sr 0.4 2~10000Li 0.06 0.2~1000Ta 0.0050.02~50000Zr 0.015 0.5~1000*可根据含量确定称样量及稀释倍数扩大测定范围上限。
2 规范性引用文件下列文件中的条款通过本方法的本部分的引用而成为本部分的条款:下列不注日期的引用文件,其最新版本适用于本方法。
GB/T20001.4 标准编写规则第4部分:化学分析方法。
GB/T14505 岩石和矿石化学分析方法总则及一般规定。
GB6379 测试方法的精密度通过实验室间试验确定标准测试方法的重复性和再现性。
GB/T14496-93 地球化学勘查术语。
3 方法提要试料用氢氟酸、硝酸、高氯酸分解并赶尽高氯酸,用王水溶解后,移至聚乙烯试管中,定容,摇匀。
分取部分澄清溶液,用硝酸(3+97)稀释至1000倍(指试料总稀释系数为1000)后,在等离子体质谱仪上测定。
4 试剂除非另有说明,在分析中仅使用确认为优级纯的试剂和去离子水。
在空白试验(6.2)中,若已检测到所用优级纯试剂中含有大于以上元素方法检出限的含量,并确认已经影响试料中以上元素低量的测定,应净化试剂。
4.1 盐酸(1.19g/mL)4.2 硝酸(1.40g/mL)4.3 氢氟酸(1.13g/mL)4.4 高氯酸(1.67g/mL)4.5 硝酸(3+97)4.6 王水750mL盐酸(4.1)与250mL硝酸(4.2)混合,摇匀。
- 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。
9种谱校正方法及matlab 程序代码采样间隔归一化成1T ∆=,采样长度为N .这样FFT 离散谱线为0,1)i X i N =-(,相应的频率分辨率2/(1/)N f N ωπ∆=∆=. 设FFT 离散谱线局部极高谱线为m (为了数学上简洁,假定从0开始,注意在MA TLAB 环境下数组实际操作的是从1开始),记频偏量δωδω=∆. 我们需要使用谱线m 和与之相邻一条次高谱线,记这连续两条谱线中左边一条序号为M (当次高谱线在m 左侧时1M m =-,反之M m =).下面列出若干算法的δ计算公式1. 加矩形窗的精确谱校正[1]i i i X U jV =+111()sin()()cos()M M M M opt M MV V M U U M K U U ωω+++-∆+-∆=- 1211cos()sin()cos()sin()opt M M opt M M K M Z V U M K M Z V U M ωωωωωω++-∆⎡⎤=+⎢⎥∆⎣⎦-∆+∆⎡⎤=+⎢⎥∆+∆⎣⎦ 2121cos()cos()()Z M Z M M m Z Z ωωωδ∆+∆-∆=+-- 2. 加矩形窗情形,采用解析单频模型的幅值比校正[1, 2]11||()||||M M M X M m X X δ++=+-+ 3. 加汉宁窗情形,采用解析单频模型的幅值比校正[1, 2]112||||()||||M M M M X X M m X X δ++-=+-+ 4. 加矩形窗情形,采用解析单频模型的复比值校正[3]11Re ()M M M X M m X X δ++⎛⎫=+- ⎪-⎝⎭5. 加汉宁窗情形,采用解析单频模型的复比值校正[3]112()M M M MX X M m X X δ+++=+-- 6. 加矩形窗情形,采用解析单频模型的复合复比值校正[3]11Re ()M m M M X M m X X δ++⎛⎫=+- ⎪-⎝⎭ 11m R m mX X X δ++=-,1111m m L m m m m X X X X X X δ---=-=-- 0.5)0.5)m L m R δδδδδ=-++((7. 加汉宁窗情形,采用解析单频模型的复合复比值校正[3]112Re ()M M m M M X X M m X X δ++⎛⎫+=+-⎪-⎝⎭ 112m m R m mX X X X δ+++=-,1111221m m m m L m m m m X X X X X X X X δ----++=-=-- 0.5)0.5)m L m R δδδδδ=-++((8. 加矩形窗,Quin 校正[4]11Re()Re(),Re()Re()m m L R m m X X X X αα-+== 11L R L R L Rααδδαα==---,, 00 R R L R δδδδδ>>⎧=⎨⎩当且其它9. 加汉宁窗,Quin 校正[4]11Re()Re(),Re()Re()m m L R m m X X X X αα-+== 212111L RL R L Rααδδαα++==---,,00 R R L R δδδδδ>>⎧=⎨⎩当且其它References1. Schoukens, J., R. Pintelon,H. Van Hamme. The interpolated fast Fourier transform: Acomparative study . IEEE Transactions on Instrumentation and Measurement. 1992, 41(2):226-232.2. 谢明,丁康. 频谱分析的校正方法.振动工程学报. 1994, 7(2): 172-179.3. 陈奎孚, 王建立,张森文. 频谱校正的复比值法.振动工程学报(已投). 2007.4. Quinn, B.G. Estimating frequency by interpolation using Fourier coefficients.IEEETransactions on Signal Processing. 1994, 42(5): 1264-1268.%========================这是调用调试==================DT=1;N=1024;PHI=pi/3;Ampl=1;CiR=11.9; %cycles in recordFreq=CiR/(DT*N); %frequencyTV=[0:N-1];DatVec=Ampl*cos(Freq*TV*2*pi+PHI);FV=fft(DatVec);figuresubplot(2,1,1);plot(TV,DatV ec);subplot(2,1,2);plot(abs(FV(1:round(N/2.56))));grid on[MV,MI]=max(abs(FV));%加矩形窗的解析校正--1FreqShift=SpecCorr(FV,MI,N,1);%加矩形窗的解析单频模型校正--2FreqShift=SpecCorr(FV,MI,N,2);%加汉宁窗的解析单频模型校正--3HanDat=DatVec.*hanning(N,'periodic')';FV=fft(HanDat);FreqShift=SpecCorr(FV,MI,N,3);%加矩形窗的解析单频模型校正+复比值法--4FV=fft(DatVec);FreqShift=SpecCorr(FV,MI,N,4);%加汉宁窗的解析单频模型校正+复比值法--5HanDat=DatVec.*hanning(N,'periodic')';FV=fft(HanDat);FreqShift=SpecCorr(FV,MI,N,5);%加矩形窗的解析单频模型校正+复比值法+左右平均--6FV=fft(DatVec);FreqShift=SpecCorr(FV,MI,N,6);%加汉宁窗的解析单频模型校正+复比值法+左右平均--7HanDat=DatVec.*hanning(N,'periodic')';FV=fft(HanDat);FreqShift=SpecCorr(FV,MI,N,7);%加矩形窗的解析单频模型校正+Quinn算法--8FV=fft(DatVec);FreqShift=SpecCorr(FV,MI,N,8);%加汉宁窗的解析单频模型校正+Quinn算法--9HanDat=DatVec.*hanning(N,'periodic')';FV=fft(HanDat);FreqShift=SpecCorr(FV,MI,N,9);===========这是子程序======================%spectrum correction assemble% the sampling interval is 1 s (or unitary)%Input: SpecVec--Discrte Fourier Spectrum from FFT% PeakIdx--the peak index, noting the matrix in MatLab start from 1% TL--the length (or the point number) of the FFT% method--correction method% output: PeakShift-- the corrected peak shifting from the peak in discrete% spectrumfunction PeakShift=SpecCorr(SpecVec,PeakIdx,TL,method)% picking up the second highest spectrum lineif(abs(SpecVec(PeakIdx-1))>abs(SpecVec(PeakIdx+1)))IP=[PeakIdx-1,PeakIdx];ShiftCorr=-1; %shift aligning with the PeakIdxelseIP=[PeakIdx,PeakIdx+1];ShiftCorr=0; %shift aligning with the PeakIdxendII=IP(1)-1; % noting that the index of a matrix in MATLAB starts from 1, not zero if(method==1) %an analyitic solution for rectangular windowU=real(SpecVec(IP));V=imag(SpecVec(IP));DW=2*pi/TL;KOPT=(sin(II*DW)*(V(2)-V(1))+cos(II*DW)*(U(2)-U(1)))/(U(2)-U(1));Z=V.*(KOPT-cos((IP-1)*DW))./(sin(DW*(IP-1)))+U;Tmp1=(Z(2)*cos(DW*(II+1))-Z(1)*cos(DW*II))/(Z(2)-Z(1));PeakPos=acos(Tmp1)/DW;PeakShift=PeakPos-(PeakIdx-1);elseif(method==2) %based on the analytical-single-tone model for rectangular window PeakShift=abs(SpecVec(IP(2)))/(abs(SpecVec(IP(2)))+abs(SpecVec(IP(1))));PeakShift=PeakShift+ShiftCorr; %shift aligning with the PeakIdxelseif(method==3) %based on the analytical-single-tone model for Hanning windowPeakShift=(2*abs(SpecVec(IP(2)))-abs(SpecV ec(IP(1))))/(abs(SpecV ec(IP(2)))+abs(SpecVec(IP(1 ))));PeakShift=PeakShift+ShiftCorr; %shift aligning with the PeakIdxelseif(method==4) %based on the analytical-single-tone model for rectangular window with complex correctionPeakShift=real(SpecVec(IP(2))/(SpecVec(IP(2))-SpecVec(IP(1))));PeakShift=PeakShift+ShiftCorr; %shift aligning with the PeakIdxelseif(method==5) %based on the analytical-single-tone model for Hanning window with complex correctionPeakShift=(2*SpecVec(IP(2))+SpecVec(IP(1)))/(SpecVec(IP(2))-SpecVec(IP(1)));PeakShift=real(PeakShift)+ShiftCorr; %shift aligning with the PeakIdx elseif(method==6) %based on the analytical-single-tone model for rectangular window with complex correction+averagePeakShift=real(SpecVec(IP(2))/(SpecVec(IP(2))-SpecVec(IP(1))));MaxPeakShift=PeakShift+ShiftCorr; %shift aligning with the PeakIdxLeftShift=real(SpecVec(PeakIdx)/(SpecVec(PeakIdx)-SpecVec(PeakIdx-1)))-1;RightShift=real(SpecVec(PeakIdx+1)/(SpecVec(PeakIdx+1)-SpecVec(PeakIdx)));%averagePeakShift=(0.5-MaxPeakShift)*LeftShift+(0.5+MaxPeakShift)*RightShift;elseif(method==7) %based on the analytical-single-tone model for Hanning window with complex correction+average,????PeakShift=(2*SpecVec(IP(2))+SpecVec(IP(1)))/(SpecVec(IP(2))-SpecVec(IP(1)));MaxPeakShift=real(PeakShift)+ShiftCorr; %shift aligning with the PeakIdxLeftShift=(2*SpecVec(PeakIdx)+SpecV ec(PeakIdx-1))/(SpecVec(PeakIdx)-SpecVec(PeakIdx-1))-1;RightShift=(2*SpecVec(PeakIdx+1)+SpecV ec(PeakIdx))/(SpecVec(PeakIdx+1)-SpecVec(PeakIdx ));%averagePeakShift=(0.5-MaxPeakShift)*LeftShift+(0.5+MaxPeakShift)*RightShift;elseif(method==8) %Quinn method for the rectangular windowa1 = real(SpecV ec(PeakIdx-1)/SpecVec(PeakIdx)); %lefta2 = real(SpecV ec(PeakIdx+1)/SpecVec(PeakIdx)); %rightLeftShift = a1/(1-a1);RightShift = -a2/(1-a2);if (LeftShift>0 & RightShift>0)PeakShift = RightShift;elsePeakShift = LeftShift;endelseif(method==9) %Quinn method for the Hanning windowa1 = real(SpecV ec(PeakIdx-1)/SpecVec(PeakIdx)); %lefta2 = real(SpecV ec(PeakIdx+1)/SpecVec(PeakIdx)); %rightLeftShift = (2*a1+1)/(1-a1);RightShift = -(2*a2+1)/(1-a2);if (LeftShift>0 & RightShift>0)PeakShift = RightShift;elsePeakShift = LeftShift;endendend。