能谱谱数据分解方法研究程序核数据处理课程设计

《能谱谱数据分解方法研究程序核数据处理课程设计》由会员上传分享，免费在线阅读，更多相关内容在学术论文-天天文库。

1、%本次课程设计采用的谱数据为iaea-1995文件夹下iaearfnwTSTSPEC%里面的数据。首先来看看CALIB.ASC。READ_ME.TXT中说明了这个谱数据包含的部分峰的峰位与对应能量如下：%ChannelEnergy(keV)%301122.06%1281511.00%1661661.66%2097834.84%29511173.24%32071274.54%33531332.50%运行程序，其中参数选择为：选择傅里叶变换法平滑输入3，选择高斯滤波器输入2，然后A=1,FWHM=4，对称零面积法的参数是K=2，H=3,b=1寻出来%的峰与READ_ME.TXT中说明

2、的部分峰的峰位与对应能量数据相吻合。clc;clear;[Filename,Pathname]=uigetfile('*.*','选择谱数据');fid=fopen([PathnameFilename],'r')%fid为文件指针，r表示读操作[array,count]=fscanf(fid,'%d',[1inf]);%指定格式转换后返回给矩阵array，同时返回成功的读出的数据数量count，1表示读出一个元素到一个列向量，inf表示读到文件结束返回一个与文件数据元素相同的列向量fclose(fid);%%%%下面开始能谱平滑%%%%%%%pinghuaxuanze=input

3、('请选择平滑方法:输入1选择重心法平滑输入2选择多项式最小二乘移动平滑法输入3选择傅里叶变换法输入4选择小波变换:');%************************重心法平滑****************************if(pinghuaxuanze==1)biaoji=1;fori=1:countarray_z(i)=array(i);endw=input('inputthewidthofthefilterwindow:');%w表示w点平滑公式whilemod(w,2)==0%判断输入的数是否是奇数，不是则重新输入。w=input('i

4、nputoddnumber:');endm=floor(w/2);forj=1:mfori=1:countif(i==1)array_smooth(i)=0.5*(array_z(i)+array_z(i+1));%能谱左边界做对称镜像处理elseif(i>1&&i<(count-1))array_smooth(i)=0.25*array_z(i-1)+0.5*array_z(i)+0.25*array_z(i+1);elsearray_smooth(i)=0.5*(array_z(count)+array_z(count-1));%能谱右边界做对称镜像处理endendfori=

5、1:count%将平滑好的数据放回原数组，为下一次做好准备。array_z(i)=array_smooth(i);endendfori=1:counta1(i)=array_z(i);end%***********************重心法平滑结束***************************%***********************多项式最小二乘移动平滑法*****************elseif(pinghuaxuanze==2)biaoji=2;w=input('inputthewidthofthefilterwindow:');%w为窗口宽度fwk=

6、zeros(w,1);%存储滤波器系数，产生一个1行，w列的零矩阵；当求平滑之后谱的第i点数据时，先在原始数据第i点的左、右各取m个数据点；也就是形成一个共有w=2m+1个数据点的窗口；fori=1:wk=i-ceil(w/2);%调整，计算采用k=-m:m，窗口的中心点为ceil(w/2)点fwk(i)=[1+(15/(w^2-4))*((w^2-1)/12-k^2)]/w;%2次或3次滤波器%fwkz=2.5*(w^2-7)*((w^2-1)/12-k^2)-9*((w^2-1)*(3*w^2-7)/240-k^4);%4次或5次滤波器%fwk(i)=(1+105/(4*(w

7、^2-4)*(w^2-16))*fwkz)/w;%fwk(i)=(1+(15*(3*w^2-7)/((w^2-4)*((w^2-5)^2+4)))*((w^2-1)*(3*w^2-7)/240-k^%4))/w;%箱型滤波器endarray_z=zeros(count+2*floor(w/2),1);%先将数据全部放在array_z数组中，并将边界做镜像处理，即增加2*floor(w/2)个数据fori=1:count+2*floor(w/2)if(i<=floor(w/2))a