这个程序用来寻找输入的两组数据(x，y)的峰，并计算其半高全宽和其他参数。

开始比较容易懂，后来到了程序主体部分，看不明白，看了好长时间，可能因为知识结构不全面，看不懂了。希望再此能得到大家的帮助，多谢啦。

可能主要是程序的设计思路没有明白。

程序在下面附件中，不明白的部分用红字标出

function[output] = autopeak(x,y)

% INITIALIZATION初始化

% set detection sensitivity parameters设置灵敏度参数为黄金分割1.618

% the golden ratio seems to work exceptionally fine :)

nfwhms = 1.6180;%1.6180;半高全宽

nstds = 1.6180;%1.6180;标准偏差

% put signal to columns if not so yet检查输入信号

shape = size(x);

ifshape(2) > 1

x = x';

end

shape = size(y);

ifshape(2) > 1

y = y';

end

% smooth signal, determine its baseline and detection threshold平滑信号

y = smooth(y);

threshold = nstds * std(y);%设定阈值

% plot the signal

%画图

plot(x,y);

holdon;

%PEAKDETECTION RUNDOWN开始寻峰

% hunt the peaks down one by one and measure them依次寻峰并测量，之后从信号中删除，知道小于阈值

% then eliminate them from signal and repeat until threshold

fori = 1 : length(x)

% determine peak by global maximum找出最大值，并计算基准1(平均值)。

[peaki, ipeak] = max(y);

base(1) = mean(y);

% check if big enough, else stop

ifpeaki < threshold

break

end

% save the peak and plot it保存峰值并画出

peak(i) = peaki;

loc(i) = x(ipeak);

plot(loc(i),peak(i),'ro');

% determine fwhm

% by shifting signal twice and locating minimas on both sides

% do it in two steps: using global baseline (j=1) and fwhm-based local baseline (j=2)

%通过移动信号2次，来定位两侧的最小值：第一步，使用基准1；第二部使用FWHM的基准2.

forj = 1 : 2

% compute shifted signal by peak-baseline points

y_shift = abs(y - (peak(i) + base(j))/2);

x_shift = abs(x - loc(i));

y_shift2 = y_shift + x_shift * peak(i)/max(x_shift);

% determine both fwhm points on shifted signal by sequential

% elimination

[y_temp, ifwhms(1)] = min(y_shift2);

y_shift2(ifwhms(1)-20:ifwhms(1)+20) = y_shift2(ifwhms(1)-20:ifwhms(1)+20) + peak(i);

[y_temp, ifwhms(2)] = min(y_shift2);

% compute index fwhm difference and fwhm

idfwhms = abs(diff(ifwhms));

indfwhms = floor(nfwhms*idfwhms);

fwhm(i) = abs(x(ifwhms(2)) - x(ifwhms(1)));

% compute fwhm-based guess for local baseline

base(j+1) = (y(ipeak-indfwhms) + y(ipeak+indfwhms))/2;

end

% recompute peak value with local baseline

peak(i) = y(ipeak) - base(j);

% calculate total power in peak with local baseline

y_temp = y - base(j);

peakint(i) = sum(y_temp(ipeak-idfwhms:ipeak+idfwhms)) * diff(x(1:2));

% plot the area under the peak

fork = 0:idfwhms

plot([x(ipeak-k),x(ipeak-k)],[base(j),y(ipeak-k)],'r');

plot([x(ipeak+k),x(ipeak+k)],[base(j),y(ipeak+k)],'r');

end

% eliminate the peak from signal for rundown to proceed

y_temp = zeros(length(y),1);

y_temp(ipeak-indfwhms:ipeak+indfwhms) = y(ipeak-indfwhms:ipeak+indfwhms);

y = y - y_temp;

end

holdoff;

% OUTPUT PREPARATION

% create nice output matrix, sort peaks by location order and number them

output(:,2) = peak';

output(:,3) = loc';

output(:,4) = fwhm';

output(:,5) = peakint';

output = sortrows(output,3);

output(:,1) = (1:length(peak))';

[本帖最后由 edifiers2008 于 2009-1-21 09:49 编辑]