核数据处理课程设计能谱谱数据分解方法研究程序.docx
- 文档编号:26103069
- 上传时间:2023-06-17
- 格式:DOCX
- 页数:18
- 大小:25.23KB
核数据处理课程设计能谱谱数据分解方法研究程序.docx
《核数据处理课程设计能谱谱数据分解方法研究程序.docx》由会员分享,可在线阅读,更多相关《核数据处理课程设计能谱谱数据分解方法研究程序.docx(18页珍藏版)》请在冰豆网上搜索。
核数据处理课程设计能谱谱数据分解方法研究程序
%本次课程设计采用的谱数据为iaea-1995文件夹下iaearfnwﻩTSTSPEC
%里面的数据。
首先来看看CALIB.ASC。
READ_ME.TXT中说明了这个谱数据包含的部分峰的峰位与对应能量如下:
% Channel Energy (keV)
%ﻩ301 122.06
% 1281 511.00
%ﻩ1661 661.66
% 2097 834.84
%ﻩ 2951 1173.24
% 3207 1274.54
%ﻩ3353 1332.50
%运行程序,其中参数选择为:
选择傅里叶变换法平滑输入3,选择高斯滤波器输入2,然后A=1,FWHM=4,对称零面积法的参数是K=2,H=3,b=1寻出来
%的峰与READ_ME.TXT中说明的部分峰的峰位与对应能量数据相吻合。
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('请选择平滑方法:
\n输入1选择重心法平滑\n输入2选择多项式最小二乘移动平滑法\n输入3选择傅里叶变换法\n输入4选择小波变换:
\n');
%************************重心法平滑****************************
if(pinghuaxuanze==1)
biaoji=1;
for i=1:
count
array_z(i)=array(i);
end
w=input('input thewidth ofthefilterwindow:
');%w表示w点平滑公式
whilemod(w,2)==0 %判断输入的数是否是奇数,不是则重新输入。
w=input('input oddnumber:
');
end
m=floor(w/2);
forj=1:
m
fori=1:
count
if(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);
else
array_smooth(i)=0.5*(array_z(count)+array_z(count-1)); %能谱右边界做对称镜像处理
end
end
for i=1:
count %将平滑好的数据放回原数组,为下一次做好准备。
array_z(i)=array_smooth(i);
end
end
fori=1:
count
a1(i)=array_z(i);
end
%***********************重心法平滑结束***************************
%***********************多项式最小二乘移动平滑法*****************
elseif(pinghuaxuanze==2)
biaoji=2;
w=input('input the width ofthefilterwindow:
'); %w为窗口宽度
fwk =zeros(w,1);%存储滤波器系数,产生一个1行,w列的零矩阵;当求平滑之后谱的第i点数据时,先在原始数据第i点的左、右各取m个数据点;也就是形成一个共有w=2m+1个数据点的窗口;
fori=1:
w
k=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^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; %箱型滤波器
end
array_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))
array_z(i)=array(-i+ceil(w/2));
elseif(i>count+floor(w/2))
array_z(i)=array(-(i-floor(w/2))+2*count+1);
else
array_z(i)=array(i-floor(w/2));
end
end
a1=zeros(1,count);
for i=1:
count%做矩阵乘法,i=1:
count即完成能谱滤波
forj=1:
w
SMZ(j)=array_z(i+j-1);
end
a1(i)=SMZ*fwk;
end
%*************************多项式最小二乘移动平滑结束****************
%*******************************傅里叶变换法开始********************
elseif(pinghuaxuanze==3)
biaoji=3;
y=fft(array,count); %对原始谱进行傅里叶变换
lvboqixuanze=input('请选择滤波器:
\n输入1选择理想低通滤波器\n输入2选择高斯型滤波器:
\n');
%*************************理想低通滤波器傅立叶变换程序开始*************************
pyy=y.*conj(y); %计算y的模平方,conj为求复数的共轭
ss=0;flag=0;%flag为噪声幅度最大值ss的标志
for i=floor(count/2*3/4):
floor(count/2)
if(ss flag=i; %flag为噪声幅度最大值ss的标志 ss=pyy(i);%找出最大幅度E对应的频率w1 end end R=5.0;flag1=0; if(flag1==0) R=R-0.01; tt=R*ss; fori=floor(count/2*3/4)-1: -1: 1 if(pyy(i)>tt) flag1=i;break;%从最大频率的3/4处向前找出R*w1的点,此点为信号起主要作用的点 end end end for i=flag1: floor(count/2*3/4) %找到信号起主要作用的点,从该点向后找出第一个频率低于w1的点,此点即为切断频率MFC if(pyy(i)<ss) flag=i;break%找到信号的切断频率并放在flag中 end end if(lvboqixuanze==1) yli=y; yli(flag+1: count-flag+1)=0;%在截断频率后赋值为0,这里考虑到对称的原因 a1=real(ifft(yli,count));%反变换后求实部 %***************************理想低通滤波器傅立叶变换程序结束************************ %*************************高斯型滤波器开始************** else(lvboqixuanze==2) A=input('请输入A: \n'); FWHM=input('请输入半高宽: \n'); delta=FWHM/2.355; %delta为高斯宽度 fori=1: ceil(count/2)%ceil表示取整 if(i>flag) ygao(i)=0;ygao(count-i)=0; else w=2*pi*i/count; ygao(i)=y(i)*A*exp(-0.5*delta^2*w^2); ygao(count-i+1)=y(count-i+1)*A*exp(-0.5*delta^2*w^2); end end gd=real(ifft(ygao,count));%反变换后求实部 fori=1: count a1(i)=gd(i); end end %************************傅立叶变换法结束**************** %***********************小波变换开始******************* else(pinghuaxuanze==4) biaoji=4; [c,l]=wavedec(array,3,'sym8'); a3=appcoef(c,l,'sym8',3); %提取小波低频系数 d3=detcoef(c,l,3); %提取小波高频系数 d2=detcoef(c,l,2); d1=detcoef(c,l,1); delta=median(abs(c))/0.6745; THR=delta*sqrt(2*log(count)); hardd1=wthresh(d1,'h',THR); %均为硬阈值,大于阈值的加'h',小于阈值的减‘h' hardd2=wthresh(d2,'h',THR); hardd3=wthresh(d3,'h',THR); c2=[a3hardd3hardd2hardd1]; a1=waverec(c2,l,'sym8');%反变换回来 %*************************小波变换结束***************** end %***************平滑结束****************************** plot(array,'r-*'); xlabel('道址'); ylabel('能量'); holdon; ifbiaoji==1 plot(a1,'k-'); elseifbiaoji==2 plot(a1,'k-'); elseif biaoji==3 plot(a1,'k-'); elseifbiaoji==4 plot(a1,'k-'); end title('平滑前与平滑后的谱'); text(4000,8000,'\leftarrow平滑前后对比','FontSize',20); holdoff; fori=1: count; ifa1(i)==0; a1(i)=1; end end %*****************下面开始峰位确定********************* %%%%%%%%%%%%%%%%%%%先采用对称零面积寻峰法(矩形波函数)**** disp('下面开始输入对称零面积法寻峰'); disp('下面开始输入对称零面积法的各参数'); disp('如果是方波的话有k=1'); K=input('请输入参数k=? : \n'); H=input('请输入参数半宽度H=? (正奇数): \n'); m=((2*K+1)*H-1)/2; w=2*m+1; b=input('请输入参数b=? : \n'); a=2*K*b; %K=4; %H=2*K+1; %w=3*H; %b=1; %a=2*K*b; m1=floor(w/2); temporary=zeros((count+2*m1),1); %以下循环为镜像处理数据 fori=1: count+2*m1 if(i<=m1); temporary(i)=a1(ceil(w/2)-i); elseif(i>(count+m1)) temporary(i)=a1(-(i-m1)+2*count+1); else temporary(i)=a1(i-m1); end end A=zeros(count,1); %以下循环为对称零面积褶积变换 fori=ceil(w/2): count+m1; forj=-(w-1)/2: (w-1)/2; if abs(j)<=(H-1)/2; T=a; else T=-b; end A(i-m1,1)=A(i-m1,1)+T*temporary(i+j); end end fori=1: count; %数据转制 SSiFENZI(i,1)=A(i,1); end B=zeros(count,1); fori=ceil(w/2): count+m1; forj=-(w-1)/2: (w-1)/2; ifabs(j)<=(H-1)/2; T=a^2; else T=b^2; end B(i-m1,1)=B(i-m1,1)+T*temporary(i+j); end end fori=1: count; %数据转存 SSiFENMU(i,1)=B(i,1); end fori=1: count; %计算SS SS(i,1)=SSiFENZI(i,1)/sqrt(SSiFENMU(i,1)); end p=1; q=1; f=30; %灵敏因子或称为找峰阈值 fori=1: count; %寻峰 ifSSiFENZI(i)<0; fpdatablow(p,1)=i; fpdatablow(p,2)=SSiFENZI(i); p=p+1; elseifSS(i)>f; fpdataup(q,1)=i; fpdataup(q,2)=SSiFENZI(i); q=q+1; end end p=1; fori=2: length(fpdataup(: 1))-1; %确定峰位 iffpdataup(i,2)>fpdataup(i+1,2)&&fpdataup(i,2)>fpdataup(i-1,2); mpeak(p,1)=fpdataup(i,1); p=p+1; end end for i=1: length(mpeak(: 1));%确定边界道 j=mpeak(i); t=mpeak(i); peak(i)=t+(a1(t+1)-a1(t-1))/(2*a1(t)-a1(t+1)-a1(t-1))/2; %*******************************确定左右边界道****************** %*********************************采用0.2倍峰高法确定边界******** whilea1(j)>0.2*a1(mpeak(i)); j=j-1; end pboard(i,1)=j; %左边界道 zuobianjie(i)=pboard(i,1); j=mpeak(i); while a1(j)>0.2*a1(mpeak(i)); j=j+1; end pboard(i,2)=j;%右边界道 youbianjie(i)=pboard(i,2); end for i=1: length(mpeak) fk(i)=pboard(i,2)-pboard(i,1); end disp('各峰的道址如下: ') sprintf('%d',mpeak) disp('峰所在道址对应的计数: ') fori=1: length(mpeak) sprintf('%d',a1(mpeak(i))) end disp('各峰的准确道址如下: ') sprintf('%d',peak) disp('各峰的左边界道址如下: ') sprintf('%d ',zuobianjie) disp('各峰的右边界道址如下: ') sprintf('%d',youbianjie) %**********下面开始能量刻度************* disp('如果你是用高斯滤波器A=1,FWHM=4,对称零面积法的参数是K=2,H=3,b=1的话,寻得的峰恰好READ_ME.TXT所示(这里是特殊情况)') biaoji3=input('各参数是否如上所述? 是的话请输1,不是的话请输2\n'); ifbiaoji3==1; nengliang=[122.06511.00661.66834.841173.241274.541332.50]; p=polyfit(peak,nengliang,1); disp('二次多项式y=kx+b的拟合系数为: 这里第一个数据是k,第二个数据是b,x对应道址,y对应能量'); sprintf('%f',p) else daozhi=input('请输入寻得的峰的道址若干: 用[]形式输入数据\n'); nengliang=input('请输入道址对应的能量: 用[]形式输入数据单位kev\n'); p=polyfit(daozhi,nengliang,1); disp('二次多项式y=kx+b的拟合系数为: 这里第一个数据是k,第二个数据是b,x对应道址,y对应能量'); sprintf('%f ',p) end fori=1: count; y(i)=p (1)*i+p(2); end figure; plot(y); xlabel('道址'); ylabel('能量kev'); title('刻度好的能量刻度方程曲线'); %*****************以下是画寻峰之后的得图*************** figure; pp=plot(a1); set(pp,'Color','blue','LineWidth',3); xlabel('道址'); ylabel('计数'); title('寻峰之后的道址与计数的关系图'); holdon; hh=axis; for(i=1: length(peak)) plot([peak(i),peak(i)],[hh(3),hh(4)],'r'); end holdoff; %*************************************************** biaoji4=input('下面计算峰面积,输入1采用瓦森峰面积法,输入2采用总峰面积法'); if biaoji4==1; %******************下面计算峰面积(这里采用瓦森峰面积法)******************** disp('计算峰面积: 下面采用瓦森峰面积法'); n=2;%经调试,取2时deta最小 y=zeros(1,10000); s=zeros(1,length(mpeak)); for i=1: length(mpeak) B1=zeros(1,length(mpeak)); B2=zeros(1,length(mpeak)); if(fk(i)/2>n) B1(i)=(a1(youbianjie(i))-a1(zuobianjie(i)))*(mpeak(i)-zuobianjie(i)-n)/fk(i)+a1(zuobianjie(i)); B2(i)=(a1(youbianjie(i))-a1(zuobianjie(i)))*(mpeak(i)-zuobianjie(i)+n)/fk(i)+a1(zuobianjie(i)); forj=(mpeak(i)-n): (mpeak(i)+n) y(i)=y(i)+a1(j); end else disp('您输入的n值不满足fk(i)/2>n的条件'); break; end s(i)=y(i)-(2*n+1)*(B1(i)+B2(i))/2; end elsebiaoji4==2; %******************************************下面用总峰面积法********************* disp('计算峰面积: 下面采用总峰面积法'); y1=zeros(1,length(mpeak)); fori=1: length(mpeak) forj=(zuobianjie
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 数据处理 课程设计 能谱谱 数据 分解 方法 研究 程序