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

％本次课程设计采用的谱数据为ｉ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:m

fｏ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);

else

arraｙ_smootｈ（i）=0．5*(array_z(coｕnt)+array_z(ｃounｔ－１)); %能谱右边界做对称镜像处理

end

ｅnd

foｒi=1:count %将平滑好的数据放回原数组，为下一次做好准备。

arraｙ＿ｚ(ｉ)＝ａrray＿smootｈ（i);

end

end

ｆor i＝1:cｏunt

a1(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：w

k=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);