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

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

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

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

4、('input odd number:');endm=floor(w/2);for j=1:mfor i=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);else% 能谱右边界做对称镜array_smooth(i)=0.5*(array_z(count)+array_z(coun

5、t-1);像处理end endfor i=1:count%将平滑好的数据放回原数组，为下一次做好准备。array_z(i)=array_smooth(i); endendfor i=1:count a1(i)=array_z(i); end%* 重心法平滑结束*%* 多项式最小二乘移动平滑法*elseif(pinghuaxuanze=2) biaoji=2; w=input('input the width of the filter window:');%w 为窗口宽度fwk = zeros(w,1); %存储滤波器系数，产生一个1 行， w 列的零矩阵；当求平滑之后谱的第i

6、 点数据时，先在原始数据第i 点的左、右各取m 个数据点；也就是形成一个共有 w=2m+1 个数据点的窗口； for i=1:w k=i-ceil(w/2); % 调整，计算采用 k=- m:m ，窗口的中心点为 ceil( w/2) 点 fwk(i)=1+(15/(wA2-4)*(wA2-1)/12-kA2)/w;%2 次或 3 次滤波器%fwkz=2.5*(wA2-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+(

7、15*(3*wA2-7)/(wA2-4)*(wA2-5)A2+4)*(wA2-1)*(3*wA2-7)/240-kA%4)/w;%箱型滤波器endarray_z = zeros(count+2*floor(w/2),1); %先将数据全部放在 array_z 数组中，并将边界做镜像处理，即增加2*floor(w/2) 个数据for i=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*

8、count+1); elsearray_z(i)=array(i-floor(w/2); end enda1=zeros(1,count);for i=1:count%做矩阵乘法 ,i=1:count 即完成能谱滤波for j=1:wSMZ(j)=array_z(i+j-1);enda1(i)=SMZ*fwk; end%*%*elseif(pinghuaxuanze=3)biaoji=3;y=fft(array,count);多项式最小二乘移动平滑结束傅里叶变换法开始%对原始谱进行傅里叶变换*lvboqixuanze=input(' 请选择滤波器 :n 输入 1 选择理想低通滤波器n

9、输入 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<pyy(i)flag=i; %flag 为噪声幅度最大值ss 的标志ss=pyy(i);% 找出最大幅度E 对应的频率w1endendR=5.0;flag1=0;if(flag1=0)R=R-0.01;tt=R*ss;for i=floor(coun

10、t/2*3/4)-1:-1:1 if(pyy(i)>tt)flag1=i;break; %从最大频率的 3/4 处向前找出 R*w1 的点，此点为信号起主要作用的点endendendfor i=flag1:floor(count/2*3/4)%找到信号起主要作用的点，从该点向后找出第一个频率低于w1 的点，此点即为切断频率MFCif(pyy(i)<ss)flag=i;break%找到信号的切断频率并放在flag 中endendif (lvboqixuanze=1)yli=y;0，这里考虑到对称的原因yli(flag+1:count-flag+1)=0;% 在截断频率后赋值为a1=r

11、eal(ifft(yli,count);% 反变换后求实部%*理想低通滤波器傅立叶变换程序结束*高斯型滤波器开始*else (lvboqixuanze=2) A=input(' 请输入 A:n'); FWHM=input(' 请输入半高宽 :n');delta=FWHM/2.355; %delta 为高斯宽度 for i=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*deltaA2*wA2

12、); ygao(count-i+1)=y(count-i+1)*A*exp(-0.5*deltaA2*wA2); end end gd=real(ifft(ygao,count);% 反变换后求实部 for i=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=det

13、coef(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=a3 hardd3 hardd2 hardd1; a1=waverec(c2,l,'sym8');% 反变换回

14、来 %*小波变换结束*end %* 平滑结束 * plot(array,'r-*'); xlabel(' 道址 '); ylabel(' 能量 '); hold on; if biaoji=1 plot(a1,'k-'); elseif biaoji=2plot(a1,'k-');elseif biaoji=3plot(a1,'k-');elseif biaoji=4plot(a1,'k-');endtitle(' 平滑前与平滑后的谱');text(4000,8000

15、,'leftarrow 平滑前后对比 ','FontSize',20);hold off;for i=1:count;if a1(i)=0;a1(i)=1;endend if abs(j)<=(H-1)/2;%*卜面开始峰位确定*% 先采用对称零面积寻峰法（矩形波函数）disp（'下面开始输入对称零面积法寻峰）；disp（'下面开始输入对称零面积法的各参数）;disp（'如果是方波的话有 k=1'）;K=input（' 请输入参数k=?:n'）;H=input（' 请输入参数半宽度H=? （正奇数）

16、: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);for i=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);elsetemporary(i)=

17、a1(i-m1);endendA=zeros(count,1);for i=ceil(w/2):count+m1;for j=-(w-1)/2:(w-1)/2;% 以下循环为对称零面积褶积变换T=a;elseT=-b;endA(i-m1,1)= A(i-m1,1)+T*temporary(i+j);endendfor i=1:count;%数据转制SSiFENZI(i,1)=A(i,1);endB=zeros(count,1);for i=ceil(w/2):count+m1;for j=-(w-1)/2:(w-1)/2;if abs(j)<=(H-1)/2;T=aA2;elseT=bA

18、2;endB(i-m1,1)=B(i-m1,1)+T*temporary(i+j);endendfor i=1:count;%数据转存SSiFENMU(i,1)=B(i,1);endfor i=1:count;%计算 SSSS(i,1)=SSiFENZI(i,1)/sqrt(SSiFENMU(i,1);endp=1;q=1;f=30;%灵敏因子或称为找峰阈值for i=1:count;%寻峰if SSiFENZI(i)<0;fpdatablow(p,1)=i;fpdatablow(p,2)=SSiFENZI(i);p=p+1;elseif SS(i)>f;fpdataup(q,1)

19、=i;fpdataup(q,2)=SSiFENZI(i);q=q+1;end endp=1;%确定峰位for i=2:length(fpdataup(:,1)-1;if fpdataup(i,2)>fpdataup(i+1,2)&&fpdataup(i,2)>fpdataup(i-1,2);mpeak(p,1)=fpdataup(i,1);p=p+1;endendfor 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(

20、t-1)/2;%*%*确定左右边界道*采用 0.2 倍峰高法确定边界*while a1(j)>0.2*a1(mpeak(i);j=j-1;endpboard(i,1)=j;%左边界道zuobianjie(i)=pboard(i,1);j=mpeak(i);while a1(j)>0.2*a1(mpeak(i);j=j+1;endpboard(i,2)=j; %右边界道youbianjie(i)=pboard(i,2);end for i=1:length(mpeak)fk(i)=pboard(i,2)-pboard(i,1);enddisp('各峰的道址如下：)sprint

21、f('%d ',mpeak)disp('峰所在道址对应的计数：')for i=1:length(mpeak)sprintf('%d ',a1(mpeak(i)enddisp('各峰的准确道址如下：)sprintf('%d ',peak)disp('各峰的左边界道址如下：)sprintf('%d ',zuobianjie)disp('各峰的右边界道址如下：) sprintf('%d ',youbianjie)%*disp（'如果你是用高斯滤波器A=1,FWHM=4 ,对

22、称零面积法的参数是K=2 , H=3,b=1的话,寻得的峰恰好READ_ME.TXT 所示 （这里是特殊情况）'）biaoji3=input（' 各参数是否如上所述？是的话请输 1，不是的话请输2n'）;if biaoji3=1;nengliang=122.06 511.00 661.66 834.84 1173.24 1274.54 1332.50;p=polyfit(peak,nengliang,1);disp('二次多项式y=kx+b的拟合系数为：这里第一个数据是k,第二个数据是b, x对应道址，y 对应能量 ');sprintf('%f

23、',p)else daozhi=input('请输入寻得的峰的道址若干：用口形式输入数据n');nengliang=input(' 请输入道址对应的能量：用 形式输入数据单位 kevn');p=polyfit(daozhi,nengliang,1);disp('二次多项式y=kx+b的拟合系数为：这里第一个数据是k,第二个数据是b, x对应道址， y 对应能量 ');sprintf('%f ',p) end for i=1:count;y(i)=p(1)*i+p(2);end figure; plot(y); xlabel

24、(' 道址 '); ylabel(' 能量kev');title(' 刻度好的能量刻度方程曲线 ');%* 以下是画寻峰之后的得图 * figure; pp=plot(a1); set(pp,'Color','blue','LineWidth',3); xlabel(' 道址 '); ylabel(' 计数 '); title(' 寻峰之后的道址与计数的关系图 '); hold on; hh=axis; for(i=1:length(peak)plot(peak(i),peak(i),hh(3),hh(4),'r');endhold off;%*2 采用总峰面积法'）;*biaoji4=input(' 下面计算峰面积，输入 1 采用瓦森峰面积法，输入if biaoji4=1;%*