array_smooth(i)=0.2438*array(i)+0.205*(array(i+1)+array(i+1))+0.1218*(array(i+2)+array(i+2))+0.05126*(array(i+3)+array(i+3)); %对称做镜像处理
end
elseif(i>k1&i<(count-k1))
array_smooth(i)=0.2438*array(i)+0.205*(array(i-1)+array(i+1))+0.1218*(array(i-2)+array(i+2))+0.05126*(array(i-3)+array(i+3)); else
array_smooth(i)=0.2438*array(i)+0.205*(array(i-1)+array(i-1))+0.1218*(array(i-2)+array(i-2))+0.05126*(array(i-3)+array(i-3)); end end
x2=1:count;
plot(x2,array_smooth)
p=ones(1,count);
m=1; %控制参数1 j=2; %控制参数2
for i=1+j:count-j-m
p(i)=(array_smooth(i)*array_smooth(i+m))/(array_smooth(i-j)*array_smooth(i+m+j)); end
k=1.5; %峰判断控制参数 for i=1:count
upper(i)=1+k/sqrt(array_smooth(i)); lower(i)=1-k/sqrt(array_smooth(i)); end%简单比较法寻峰 factor=1; for i=1+factor:count-factor
if(p(i)>upper(i)&p(i)>p(i+1)&p(i)>p(i-1)) peak_position=[peak_position i]; end
end
plot(array_smooth,'b*'); hold on;
peak=array(peak_position) plot(peak_position,peak,'r+'); 实验结果如下:
四、本文第二种思想:gauss8阶函数拟合法,此方法的主要思想是:不将数据平滑处理自然就没有了平滑公式,直接用八阶高斯函数拟合源数据,在对拟合后的函数进行求导,一阶导,峰位处由正到负过零点,边界处由负到正过零点 ;二阶导,峰位处负的局部最小值,边界道为正的极大值;三阶导,峰位处由负到正过零点,边界处由正到负过零点。其他过程都与前面的计算方法相同。
开始 获取能谱数据 高斯函数法对能谱进行拟合 对拟合后的函数求导 导数法的寻峰 得到计算峰面积 结束
五、MATLAB实验程序如下: % y 原始数据,1025*1 % x = (1:1025)';
% y2为前600个数据 也可用全部数据 y2 = y(1:600); % 这四个峰出现在y(1:600)中
% x2 = (1:600)'; x2 = (1:length(y2))';
% Set up fittype and options. ft = fittype( 'gauss8' ); opts = fitoptions( ft ); opts.Display = 'Off';
opts.Lower = [-Inf -Inf 0 -Inf -Inf 0 -Inf -Inf 0 -Inf -Inf 0 -Inf -Inf 0 -Inf -Inf 0 -Inf -Inf 0 -Inf -Inf 0];
opts.StartPoint = [1236 282 3.92319043722374 907 336 4.22548549053934 338.232193942756 342 5.30830398585762 308 230 6.26590931040574 254.442173171231 277 6.49267774242084 135.84970058518 329 7.80555875890908 112 1 21.9443863743451 102 398 10.6865668950446];
opts.Upper = [Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf];
% Fit model to data.
[fitresult, gof] = fit( x2, y2, ft, opts );
% Plot fit with data. figure( 'Name', 'fit gauss8' ); h = plot( fitresult, x2, y2 );
legend( h, 'y2 vs. x2', 'fit gauss8', 'Location', 'NorthEast' ); % Label axes xlabel( 'x2' );
ylabel( 'y2' );
title('gauss8拟合的效果图'); grid on
dd111 = diff(fitresult(x2));
x20 = x2(dd111(1:end-1).*dd111(2:end)<0) %导数值跨过0的值
% x20 = [228; 238; 259; 267; 281; 308; 336; 376;
figure
plot(x2,fitresult(x2),x20,fitresult(x20),'o') title('边界道址、峰位');
xz11 = 228 - (238-228); % 第一个峰左侧的边界道址(对称点) Z1 = trapz(x2(xz11:238),y2(xz11:238)) % 第一个峰的面积
Z2 = trapz(x2(267:308),y2(267:308)) % 第二个峰的面积
Z3 = trapz(x2(308:376),y2(308:376)) % 第三个峰的面积
xz42 = 397 + (397-376); % 第四个峰右侧的边界道址(对称点) Z4 = trapz(x2(376:xz42),y2(376:xz42)) % 第四个峰的面积
六、实验结果输出:
397];