peak addition/ arrary addition

I have a question about the addition of peaks in the form of 2D arrays, say, each peak is a 400X2 array.

Peak 1:

P1x = [1500:0.01:1504]' % 400 points centered at 1502

P1Y = 5.*exp(-(1-x./1502).^2/0.3) % some function of x

similarly, peak2:

P2x = [1501.8:0.01:1505.8]' % 400 points centered at 1503.8

P2Y = 12.*exp(-(1-x./1503.8).^2/0.5) % some function of x

Peak 3:

P3x = [1505:0.01:1509]' % 400 points centered at 1507

P3Y = 17.*exp(-(1-x./1507).^2/0.8)

.....etc. How to do the peak addition if there is any overlap and store them in a 2D array, so X of the peak sum will be [1501:0.01:1508] ( not from 1500 to 1509); Y of the sum is just the addition of the peaks.

what if I have hundreds of peaks? Thanks

Answer 1

[Iain]

If all the peaks have the same parameterised formula

 f =  @(A,B,C,x)(A.*exp(-(1-x./B).^2./C);
 for i = 1:numel(x)
  y(i) = sum(f(A,B,C,x(i)));
 end

Its not quick, but it will do the job for 100s of thousands of peaks. - You'll need to supply "x".

Alternatively, a lot quicker, if you use "bsxfun", and a bit of intelligence, you'll be able to get rid of the loop. - But you won't be able to handle quite so many peaks.

[Chong Tao]

Thanks a lot lain, the peaks are calculated from some pretty complex rountine. attached is the code. as you can see, it generates 3 peaks with different peak centerres. how can I add them up?

% clear all

p = 0.2;

gD=[0.015,0.02,0.031];

gL=[0.2*p, 0.32*p, 0.15*p];

v0 = [5326.65,5327.4,5328.56]; % Peak center

sigma1 = 0; sigma2 = 0; sigma3 = 0;

for k =1:length(v0)

vt = v0(k)-1:0.001:v0(k)+1; % this is the X axis for each peak

x= 2*sqrt(log(2)).*(vt-v0(k))/gD(k);

y = 2*sqrt(log(2))*gL(k)/gD(k);

a = 1/100;

Nn = 8000;

    for n = 1:1:Nn

        sgma1 = 1/(a^2*n^2 + y^2).* exp(-(a^2*n^2 + x.^2));

        sgma2 = 1/(a^2*n^2 + y^2).* exp(-(a*n + x).^2);

        sgma3 = 1/(a^2*n^2 + y^2).* exp(-(a*n - x).^2);

        sigma1 = sigma1 + sgma1;

        sigma2 = sigma2 + sgma2;

        sigma3 = sigma3 + sgma3;

    end

     Vgt = exp(-x.^2).*erfcx(y).*cos(2*x.*y)+...
            2*a*x.*sin(x.*y)/pi.*exp(-x.^2).*sin(x*y)/(x*y)+...
            2*a/pi.*(-y.*cos(2*x.*y).*sigma1 + y/2.*sigma2 +y/2.*sigma3);

figure(k);

plot(vt,Vgt,'.-b')

end

[Iain]

If the equation (not the numbers) for each peak is the same, or can be put into the same terms, then the approach will work. You simply need a set of arrays or vectors for each parameter (which I thought were amplitude (A), scaling (B), and sigma (C) for the "gaussian" equation you gave) If your peak is something funkier, say, sin(A*x^B+C)/(Dx^E + F), it'll work with that too.

 eg. 
 x = 1502:0.1:1506;
 for i = 1:number_of_peaks
  Peak(i,:) = A(i).* exp(-(1-x./B(i)).^2./C(i));
 end

Or vectorise it, using a bit more intelligence than I've applied, using bsxfun

Or:

 x = 1502.5;
 Total_at_x = sum(A.* exp(-(1-x./B).^2./C));

[Chong Tao]

Thanks again lain. I'll try what you suggested.

[Chong Tao]

I tried the above method. It works fine when the range of x is small, but it wont work if X range is big, say from 1000 to 1500, with a stepsize of 0.001. ieadlly, I'd like to calculate the peak with only +/-200 points arround the peak center and sum all the peaks.