％本次课程设计采用的谱数据为ｉaｅa－199５文件夹下ｉａearｆnwﻩTSＴSPEC%里面的数据。
首先来看看CＡLIB.ASC。
RＥAD_ME.TＸT中说明了这个谱数据包含的部分峰的峰位与对应能量如下：% Ｃhａnnｅl Ｅnergy(kｅV)％ﻩ３0112２.06% 1281 ５１1.00%ﻩ16６1 6６1.66% ２０9７834.8４％ﻩ2９５1 １173．24% ３207 1２74.54%ﻩ33５3 1332.50%运行程序,其中参数选择为:选择傅里叶变换法平滑输入3,选择高斯滤波器输入2，然后Ａ=1,FＷHM=4,对称零面积法的参数是K＝2,H＝3,b＝1寻出来%的峰与RＥAD＿ＭE．TXT中说明的部分峰的峰位与对应能量数据相吻合。
clｃ；cleａr；[Fiｌenamｅ，Pathｎame］=uigetfiｌｅ('*．*'，'选择谱数据');fid＝foｐｅｎ(［Pathnaｍe Fｉlenamｅ]，'r＇）%fid为文件指针，r表示读操作[ａrｒay,counｔ]=fsｃanf(fid,'％d',[1 inf])；%指定格式转换后返回给矩阵arｒay，同时返回成功的读出的数据数量count,１表示读出一个元素到一个列向量,inf表示读到文件结束返回一个与文件数据元素相同的列向量fｃlose(fid)；%％%％下面开始能谱平滑%%%%％％％pinghuaｘuａnze=iｎｐuｔ('请选择平滑方法:＼n输入１选择重心法平滑\ｎ输入2选择多项式最小二乘移动平滑法＼ｎ输入3选择傅里叶变换法\n输入4选择小波变换：\ｎ'）;%********＊*************＊*重心法平滑＊＊****＊***＊＊*****＊**＊******＊iｆ(pｉnghuaxuａnｚe＝=1)ｂｉaoji=1;for i=1:coｕnｔarray＿ｚ(i)＝arｒay(i)；ｅndｗ=ｉnｐut（'ｉnpｕt tｈe widtｈｏf the fｉlter ｗｉndoｗ:＇); %w表示w 点平滑公式while ｍoｄ（w,2)==0 %判断输入的数是否是奇数,不是则重新输入。
w=ｉnput(＇inｐut odｄnuｍbｅｒ:');eｎｄm＝flooｒ(w/２）;for j=1:mfｏr ｉ=１:cｏunｔｉｆ(i=＝1)array_sｍooth（ｉ)＝０.5*（aｒray＿z(ｉ)+arraｙ_z(i＋1)）; %能谱左边界做对称镜像处理ｅｌseif（ｉ＞1&&i<(cｏuｎｔ-1)）ａrray＿ｓmooth（ｉ)=0.２５＊arrａｙ＿z（ｉ-１)+0.5*aｒray_z(ｉ)+0．2５*arｒay_ｚ（i＋1);elsearraｙ_smootｈ（i）=0．5*(array_z(coｕnt)+array_z(ｃounｔ－１)); %能谱右边界做对称镜像处理endｅndfoｒi=1:count %将平滑好的数据放回原数组，为下一次做好准备。
arraｙ＿ｚ(ｉ)＝ａrray＿smootｈ（i);endendｆor i＝1:cｏunta1(i)=aｒray_z(i);ｅnd%＊＊＊＊***＊*＊*****＊＊**＊***重心法平滑结束*******＊***＊*＊*********＊＊＊＊%＊＊*＊***＊＊**＊***********多项式最小二乘移动平滑法******＊****＊**＊**ｅｌseif(pｉｎgｈuaxuanｚe=＝2)bｉaｏji=2；ｗ=inpuｔ('ｉnｐｕｔｔｈｅwiｄｔh oｆthe ｆｉlter winｄｏw:');％ｗ为窗口宽度fwk= zｅrｏs（w,1); %存储滤波器系数，产生一个1行，w列的零矩阵;当求平滑之后谱的第i点数据时,先在原始数据第i点的左、右各取m个数据点;也就是形成一个共有ｗ=2m+1个数据点的窗口;fｏr i=1：wk=i-ceｉl(w/２）; %调整,计算采用ｋ＝- m：m，窗口的中心点为cｅiｌ( w/2)点fwｋ（i）=[1＋(15/（w^2-4)）＊((w^2-1)/１2-k^2）］／w；%２次或3次滤波器%fwkz=2．5*（w^2-７）＊((ｗ^2-１)/1２-ｋ^2)-9*（(w＾2-1)*(３*ｗ^２-7）/240-k^4）; %4次或5次滤波器%fwｋ(i)=（1+1０５／(４*(ｗ^2-4）*(ｗ^2-16))＊fｗkz）／w;%ｆwk(i)=（1＋(１5*(3*ｗ＾2－７）/（(w＾2-4)*(（w^2-５)^２+4)））＊（(w^2-１)*(3*w＾2-7)/2４０-k^%4))/w；%箱型滤波器eｎｄarray＿z = zｅros(ｃouｎt+2*ｆlｏor(w/２),1);%先将数据全部放在aｒray＿z数组中,并将边界做镜像处理，即增加2*flooｒ(w/2)个数据for i=１:count+２*flｏoｒ（ｗ/２)if(ｉ<=fｌoor(w/2）)aｒray_z(i)=arｒａy(－i＋ceiｌ（ｗ/2)）；elｓｅiｆ(ｉ>ｃount＋floｏr(ｗ/２）)aｒray_z(i)=arｒay(-(i－ｆloｏr(w／2）)+2*count+1);elsｅarray_ｚ（i)=ａｒrａｙ(ｉ-flｏｏｒ(ｗ/2));enｄenda1=zerｏｓ(1,coｕｎt);fｏrｉ=1:couｎt ％做矩阵乘法,ｉ=１:cｏuｎt 即完成能谱滤波fｏr j=1:wSMZ(j)=array_z(ｉ+j－1);ｅｎda１(i)＝SＭZ*ｆwk;end%*＊*******＊********＊****＊*多项式最小二乘移动平滑结束**＊*＊＊＊＊**＊***** ％*＊***＊＊*******＊＊**＊***＊＊＊*＊****傅里叶变换法开始*＊**＊＊*＊*＊*******＊**eｌsｅｉf（pingｈuaｘuａｎzｅ==3)biaoji＝３；y=fft(arｒａy,count); %对原始谱进行傅里叶变换ｌvboqixｕaｎzｅ=input('请选择滤波器:\n输入1选择理想低通滤波器\n输入2选择高斯型滤波器:\n');%***＊****＊*****＊＊*＊＊**＊**＊理想低通滤波器傅立叶变换程序开始***********＊******＊**＊*＊＊pｙｙ=ｙ.＊cｏnj(ｙ)；%计算y 的模平方，conｊ为求复数的共轭ｓｓ=０;ｆlａg=0; %ｆｌag 为噪声幅度最大值ss 的标志for i＝ｆlｏoｒ(count/2＊3/4):floor（couｎt/2)if（sｓ<pyy(i)）flaｇ=i;%flａg为噪声幅度最大值ss 的标志ss＝pyｙ(ｉ)；%找出最大幅度Ｅ对应的频率w1enｄｅndR=5.0;flａg1=0；ｉf(flag１=＝０)R＝R-0.01;ｔt=R*ss;for i=floor（cｏunｔ/2*３/4)-１:-1：1if（ｐyy(i)>tt）flag1=i;breａk；%从最大频率的3/4处向前找出R*w1的点,此点为信号起主要作用的点endｅnｄendfor i=flag1:floor(counｔ/２＊３/4) %找到信号起主要作用的点,从该点向后找出第一个频率低于w1的点，此点即为切断频率MFＣif(ｐyｙ(i)<ｓs）ｆlag=i;ｂreaｋ%找到信号的切断频率并放在ｆl ａg 中endｅndｉｆ(ｌvboqiｘuａnzｅ＝=1)ｙli=y;ylｉ(flag+１:count-flaｇ＋1)=0;%在截断频率后赋值为０,这里考虑到对称的原因a１=real(ｉfft（yli,count）);%反变换后求实部%*****＊***＊***＊＊************理想低通滤波器傅立叶变换程序结束**＊***＊＊**********＊＊*＊*＊%***＊*＊**＊******＊＊*＊****＊＊高斯型滤波器开始*＊********＊***elsｅ(ｌvboqiｘuanze＝=2)A=inpｕt(＇请输入A:\n')；FWHＭ=inpｕt(＇请输入半高宽:\n＇)；dｅｌta=FWＨＭ/2.355; ％delｔa为高斯宽度for i=1：ceil(count／2)%ceil表示取整ｉｆ(ｉ>flag)ygaｏ(i)＝0;ygaｏ(counｔ－i)=0;eｌｓeｗ=2*pi*i／cｏunt;ygａo(i）=y（i)*Ａ*exp（-0．5*delta＾2*w^２);yｇao(count-i＋1)=ｙ(ｃount－i＋1)*A*exp(-０.5*dｅlta^2*w^2）;endendgd=rｅaｌ(ｉfｆt(ygａo，count));%反变换后求实部for i＝１:cｏunta1(i)=gd（i)；endend%*＊＊***＊*＊**＊****＊***＊**＊傅立叶变换法结束*********＊＊＊*＊*＊％***＊****＊**＊*＊＊*****＊**小波变换开始*********＊*＊*******elｓe（ｐinghuaxuanｚｅ＝=４)ｂiaoｊi=4;[c,ｌ]=waｖｅｄec（ａｒraｙ,３，＇sym8');a３＝aｐpcｏef(c,l,'syｍ８',3); %提取小波低频系数d3＝detcｏｅｆ(c,l,3)；%提取小波高频系数d２＝detcoef(ｃ,l，２);d1=dｅｔcoef(c,l,１);ｄelta=mediａn（abs(c))/0.67４5;TＨR＝deｌｔa＊sqrt(２*loｇ(ｃoｕnt））;haｒdd１=ｗｔhreｓh(d1,'h',THR);％均为硬阈值，大于阈值的加＇h＇,小于阈值的减‘h'ｈaｒdd2=ｗtｈresh（d2,'h',THＲ）；haｒｄｄ3=wｔhresh(d３,＇ｈ',THＲ);c２=［a3 hardd３hａrdd2 hａrdｄ1];a1=waverec(c2，l，'ｓym8');％反变换回来%******＊*******＊***＊****＊＊小波变换结束***＊***＊***＊*＊*＊*ｅnd%*＊*＊***＊*＊****＊平滑结束****＊**＊＊＊****＊*******＊*******plot(array，＇r-*');xlabel('道址');ylabel('能量');holｄoｎ;if ｂiａｏｊi=＝１pｌot(a1,'ｋ-');ｅlseｉｆbｉａoji=＝2plot(a1,＇k-');eｌｓeifｂiaｏji=＝3ｐｌot(a1,'k－＇);elsｅiｆbｉaoji==4pｌｏt（a1,'k-')；eｎｄtｉtle('平滑前与平滑后的谱＇)；teｘt(4000，8000，'\leftarrow 平滑前后对比'，'ＦontＳｉze'，20);holｄoff;for ｉ=１:ｃouｎｔ；if a１(i)==0；ａ1(i)＝1；eｎdｅｎd%＊***＊**＊＊****＊＊**下面开始峰位确定*******＊********＊**** %%％%%%%%%%％％%%%%%%%先采用对称零面积寻峰法(矩形波函数）**＊* disp('下面开始输入对称零面积法寻峰'）;diｓp('下面开始输入对称零面积法的各参数');dｉsp（'如果是方波的话有k＝1')；Ｋ=input（'请输入参数k=?:\n');H=ｉnpuｔ(＇请输入参数半宽度Ｈ=?(正奇数）：＼n'）;ｍ=((2*K+1)＊H－１)/2；w=2*m＋1;b＝iｎpuｔ('请输入参数ｂ=？:\ｎ＇);a=2＊K*b;%K=4;%H=2*Ｋ+1；%w=3*Ｈ;％ｂ=1；%a=2*K*b;m1=flｏｏｒ(ｗ／2）;ｔempoｒary＝zeros(（counｔ+2＊m1),1）; ％以下循环为镜像处理数据for i=1:cｏuｎt+2*m1if（i<=m１)；temporary(i）=a1(ｃｅil(w/2）－i)；ｅlｓｅｉf(i>(count+ｍ1))temｐｏｒaｒy(i)=a1(-(i-m1)+2*couｎt+1);ｅｌseｔemporary(ｉ)=a1（ｉ-m1);eｎdendA=ｚeroｓ(ｃoｕnt，1); %以下循环为对称零面积褶积变换ｆor i＝ceil(w/2):ｃｏｕnt+m1;for j=-(w-1)/2:(w-1）/2;iｆabs(ｊ)<＝(H－１)/2;Ｔ=a；ｅｌseT=-b;endA(i-m１，1)＝A(ｉ－ｍ1,1）＋T*tｅmporarｙ(i+ｊ);enｄｅｎdｆｏｒi=1:ｃounｔ;%数据转制SＳｉＦENZI（i，1)=A(i,1);ｅndＢ=zeroｓ(cｏｕnｔ，1);for i=ceiｌ(w／2)：ｃoｕnt+m1;ｆｏｒj=-（w-1）/2:（ｗ-1）／2;if abｓ（j)<＝（Ｈ－1）/2;T＝a＾２；ｅlsｅT=b^２;ｅndＢ(i-m1,1)=B(i-m1，１)＋T*tｅmporary（i+j);ｅndｅndfoｒｉ=1：count; ％数据转存SＳiFENMＵ(i，1)=B(ｉ,1);ｅｎdfor ｉ=1:cｏｕnt; %计算SSSS(i,1)=ＳＳiFENZI(ｉ，１)/sqrｔ(SSｉＦＥNMU(ｉ,1))；endｐ＝1；q=1；f＝30；%灵敏因子或称为找峰阈值for i=1:count; %寻峰ｉf SＳｉＦENZI（i)<0；fpdatablow(ｐ,１)=ｉ；fｐdaｔablow(p,2)＝SSiFENZI(ｉ);p＝p＋1;eｌsｅif ＳS(ｉ）>f;fpdatａup(q,１)=i;fpdataｕp(ｑ,2)=SSｉFENZI(ｉ);q=ｑ+1;endeｎｄp＝１；foｒi=2：ｌｅngth(ｆpdatａuｐ(：,１)）－1；%确定峰位if fpｄataｕp(i,2)＞fpdataup(ｉ＋1，2)&&fpdａtaｕp(i,2)>fpｄatａup(i-１，２)；ｍpeak(ｐ，1)=fpdatａup（i,1）;ｐ=p＋1;enｄendforｉ=1:leｎgｔh(mpeaｋ(:,1）); ％确定边界道j=mpeaｋ(i）;t=mpeaｋ(i）;peａk(ｉ)=t＋(a1(t+1）-a1(t-１))/(2*a1(t)-a１(ｔ+1)-ａ1(t－1))/2;%**＊*＊＊*******＊＊＊＊*****＊*＊*＊*＊**确定左右边界道**＊*****＊*＊****＊*＊%*＊******＊*＊*＊**＊＊****＊＊**＊***＊***采用０.2倍峰高法确定边界*****＊**ｗhiｌe a1(j)＞0．2*a1(ｍｐeａk（i））;j=ｊ-1;ｅｎｄpboard(i，1)=j；%左边界道ｚuoｂｉanｊie(i)=pｂｏａrd(i,1）;ｊ=mｐeak（i）;while a1(j)＞0．２*ａ１(ｍｐｅak(i));j＝j+1;ｅｎdpboａrd（ｉ,２)=j; %右边界道yｏubｉanｊie(ｉ）=pｂoard(i,2)；endｆor i=1:leｎgｔh(mｐeak)fk(i)=pbｏard(i,2）-pboard(i,1）;ｅnｄdisp（'各峰的道址如下：')sｐriｎtf('%ｄ',mpeak)ｄiｓｐ（＇峰所在道址对应的计数：'）for i=１：ｌengtｈ(mpeak)ｓｐｒiｎtf('％d ',a1(mｐeak（i）))eｎddｉsp('各峰的准确道址如下:')ｓprintf('%d ',ｐeaｋ)disp(＇各峰的左边界道址如下:'）ｓpｒintf('%d＇，zｕobiaｎjie)diｓp（'各峰的右边界道址如下：')spriｎｔf(＇%ｄ',youbiaｎjie)%＊＊**＊**＊＊*下面开始能量刻度＊**＊*******＊*dｉsp('如果你是用高斯滤波器A=1，FＷＨM=4,对称零面积法的参数是K＝2,H＝3,b＝1的话，寻得的峰恰好REＡD_ME．TXＴ所示（这里是特殊情况）＇)ｂｉaojｉ3=input('各参数是否如上所述?是的话请输1,不是的话请输2\n');if biaojｉ3==１；nenｇｌianｇ=［１２2.06 ５11.０0 661.66 834．84 1１73.24 １274.５4 133２.50]；p=polyfｉt(ｐeａk,ｎenｇliａｎg,１);disｐ(＇二次多项式y=ｋｘ+ｂ的拟合系数为:这里第一个数据是k，第二个数据是b,ｘ对应道址，y对应能量');sｐｒinｔｆ('%f ',ｐ)elｓe daｏｚhｉ=iｎｐｕｔ（'请输入寻得的峰的道址若干:用[]形式输入数据\ｎ'）;nengｌiang=input（＇请输入道址对应的能量：用[］形式输入数据单位kev\ｎ＇）;p=polyｆit(ｄaozhi,ｎeｎｇliａng，1）；ｄiｓｐ（'二次多项式y=kx+ｂ的拟合系数为:这里第一个数据是k，第二个数据是b，x 对应道址，y对应能量');ｓpriｎｔf(＇%f ',p)endfor i=1:coｕｎt;ｙ(i)=p（1）＊ｉ+ｐ(２）;ｅｎdｆigure；plot(y);xｌabel（'道址');ylａｂｅl（＇能量ｋｅv');ｔiｔle('刻度好的能量刻度方程曲线'）;％＊＊**＊****＊*******以下是画寻峰之后的得图**＊******＊*****figure;ｐp＝plot（a１）;sｅt(ｐp,'Ｃolor','blue','ＬineWiｄth'，3);xlabel(＇道址＇);yｌabｅl(＇计数＇）;title('寻峰之后的道址与计数的关系图');hｏld on;hh＝aｘis;for(i=１:ｌenｇth（peak))ploｔ([ｐeａk(i),peａk(i)],［hh(3),hh(４)],＇ｒ'）;enｄhoｌd off;%*＊＊***＊*****＊＊****＊*******＊**＊＊**＊＊＊****＊***＊******biaoｊｉ4=ｉnpuｔ（'下面计算峰面积，输入1采用瓦森峰面积法,输入2采用总峰面积法'); iｆbｉaoji4==1;%*******＊*****＊**＊*下面计算峰面积（这里采用瓦森峰面积法）＊*＊＊******＊*********disｐ（'计算峰面积:下面采用瓦森峰面积法')；ｎ=２;％经调试,取2时ｄeta最小y＝zerｏｓ（1,１０00０)；ｓ=zeros(1,lenｇtｈ(mpeaｋ));foｒi=1:lengｔｈ（ｍpeak)B1＝zeros(1,leｎｇｔｈ(mｐeａｋ）)；B2=zeｒos(1,lenｇｔh(mpｅａｋ））；iｆ(ｆk（i)／2>n)B1（ｉ)=(a1(yｏuｂiaｎjie(i）)-ａ１(ｚｕobiaｎjｉe(ｉ)))*(mｐeak(i)－zｕobianｊie(i)-ｎ)/fk(i)＋a1(zｕoｂiａｎjie(i));B２(i)=（a1(yoｕbianjie(i))－a1(zｕobianｊie(i)))*(ｍpｅak(i）-zuｏbiaｎｊie(i）+n）／fk(i)+ａ1(zuobiaｎjie(i))；for j＝(mｐeａｋ（i)－n):(mｐeak(ｉ)＋n）y(i)=y(ｉ)+a１(j)；eｎdelｓeｄiｓp（'您输入的n值不满足fk(i)／2>n的条件');breａk;eｎｄｓ(ｉ)=ｙ（i)-(2*n+１)*(Ｂ1(i)+ B２(i)）／2；endelｓｅbｉａoji4＝＝2；％******＊******＊*＊＊*****＊＊****＊****＊*＊***＊**下面用总峰面积法**************＊***＊＊＊ｄｉsｐ('计算峰面积:下面采用总峰面积法'）;y1=zeros(1,lengtｈ(mpeak));foｒi=1:lengｔh（mpeaｋ)foｒj=（ｚｕｏbiａnｊiｅ(ｉ)+１):(youbianjiｅ(ｉ）-１)y1(i)=y1（ｉ)+a１(j）;endｓ（i)=y1（i)－(youbｉａnjiｅ（i)-zｕobianjie(ｉ)-1）*(ａ1(ｚuｏｂiaｎｊｉe(ｉ))+a1(youbianｊiｅ(i））)/2;eｎdenddisp('各个峰的面积为:＇);ｓprｉntf('%d '，s)。