资源描述:
《能谱谱数据分解方法研究程序核数据处理课程设计.doc》由会员上传分享,免费在线阅读,更多相关内容在工程资料-天天文库。
1、%本次课程设计采用的谱数据为iaea-1995文件夹下iaearfnwTSTSPEC%里面的数据。首先来看看CALIB.ASCoREAD_ME.TXT屮说明了这个谱数据包含的部分峰的峰位与对应能量如下:%ChannelEnergy(keV)%301122.06%1281511.00%1661661.66%2097834.84%29511173.24%32071274.54%33531332.50%运行程序,其屮参数选择为:选择傅里叶变换法平滑输入3,选择高斯滤波器输入后A=1,FWHM=4,对称零面积法的参数是K=2,H=3
2、,b=l寻出来%的峰与READ_ME.TXT屮说明的部分峰的峰位与对应能暈数据相吻介。2,然clc;clear;[Filename,Pathname]=uigetfile(,*.*,,'i^择谱数据');fid=fopen([PathnameFilename],'r')%fid为文件指针,i■表示读操作[array,count]=fscanf(fid,,%d,,[linf]);%指定格式转换后返冋给矩阵array,同时返冋成功的读出的数据数量count,1表示读出一个元素到一个列向量,inf表示读到文件结朿返冋一个与文件数据
3、元素相同的列向量fclose(fid);%%%%>'ifii开始能谱平zH-%%%%%%%pinghuaxuanze=i叩ut(*请选择平滑方法:输入1选择重心法平滑输入2选择多项式最小二乘移动平滑法输入3选择傅里叶变换法输入4选择小波变换:W);1«£••土•士•士«£••土•士■一
4、、■丿•I•I/%JfIa•士•土•士•士...j1vy.1••w1A•~if(pinghuaxuanze==1)biaoji=l;fori=l:countarray_z(i)=array(i);endw=i叩ut「i叩ut
5、thewidthofthefilterwindow:');%w表示w点平滑公式whilemod(w,2)==0%判断输入的数是否是奇数,不是则重新输入。w=input(4nputoddnumber/);endm=floor(w/2);forj=l:mfori=l:countif(i==l)array_smooth(i)=0.5*(array_z(i)+array_z(i+1));%能谱左边界做对称镜像处理elseif(i>1&&i<(count-l))array_smooth(i)=0.25*array_z(i-l)+0.5
6、*array_z(i)+0.25*array_z(i+1);elsearray_smooth(i)=0.5*(array_z(count)+array_z(count-1));%能谱右边界做对称镜像处理endendfori=l:count%将平滑好的数据放冋原数组,为下一次做好准备。array_z(i)=array_smooth(i);endendfori=l:countal(i)=array_z(i);end逐***********************重心法平滑结束***************************
7、轻***********************多项式最小二乘移动平滑法*****************elseif(pinghuaxuanze==2)biaoji=2;w=input(*i叩utthewidthofthefilterwindow:');%w为窗口宽度fwk=zeros(w,l);%存储滤波器系数,产生一个1行,w列的零矩阵;当求平滑之后谱的笫i点数据时,先在原始数据第i点的左、右各取m个数据点;也就是形成一个共有w=2m+l个数据点的窗口;fori=l:wk=i-ceil(w/2);%调整,计算采用k=-
8、m:m,窗口的中心点为ceil(w/2)点fwk(i)=[1+(15/(wA2-4))*((wA2-1)/12-kA2)]/w;%2次或3次滤波器%fwkz=2.5*(vvA2-7)*((wA2-1)/12-kA2)-9*((wA2-1)*(3*wA2-7)/240-kA4);%4次或5次滤波器%fwk(i)=(1+105/(4*(wA2-4)*(wA2-16))*fwkz)/w;%fwk(i)=(1+(15*(3*wA2-7)/((wA2-4)*((wA2-5)A2+4)))*((wA2-1)*(3*wA2-7)/240-
9、kA%4))/w;%箱型滤波器endarray_z=zeros(count+2*floor(w/2),1);%先将数据全部放在array_z数组中,并将边界做镜像处理,即增加2*floor(w/2)个数据fori=l:count+2*floor(w/2)if(i<=floor(w/2))a