How to calculate spectral overlap integral?

Share icon

Thu, 08/20/2015 - 12:30 pm

Hi,

I have two waves with normalization. I would like to calculate spectral overlap integral based on the equation in the image.

Is there a snippet I can learn from?

The equation requires more information than just the area of overlap. In that case, I just need to know how to calculate the area of the overlap.

I will modify the snippet to fit the equation later.


Thanks!
You have two integrals in the figure but we can generalize the answer by considering an integral of the function g(x) where g(x) may consist of one or more factors. For example, you can have:
g(x)=F(x)*Epsilon(x)*(x^4)


In general, you compute this type of integral using Integrate1D(). Note that in the event that one of your factors above corresponds to data in a wave, you may want to write a function that returns an interpolated value. For example, if F(x) is derived from a 1D waveF, and if waveF has appropriate set wave scaling, you can use waveF(x). If your graph is plotted as waveF vs waveX you would need to use
interp(x,waveX,waveF)


I hope this helps,

A.G.
WaveMetrics, Inc.

August 20, 2015 at 03:15 pm - Permalink

Thank you for the input. I have to admit that I cannot figure out how to do it still.

These are the parameters I have
Wave 1: 400 point
Wave 2: same as above
Both have spectral overlap.
Scale to 400 to 800 as wavelength, so point 1 = 400, point 2 = 401 and etc....
Then calculate the overlap integral J

For example, at point 1, lambda = 400nm, wave 1 = 0.3, wave 2 = 0.8
J1 at point 1 = (400nm)^4 * 0.3 * 0.8 / V_sum of wave 2
J2 at point 2 = (401nm)^4 * 0.4 * 0.9 / V_sum of wave 2
....
After that, I then add all the points from J1 to J399 to get integral J.

I have figured it out using MS exce, but I don't know how to do it in IGOR.
I have tried this:

Function SpectralIntegral(w1,w2,wave_nm)
wave w1,w2,wave_nm
wavestats /q w2
variable xx = v_sum

variable j1, j2, j3,J_total
j1 = w1[0] *w2[0] * wave_nm[0]^4/xx
j2 = w1[1] *w2[1] * wave_nm[1]^4/xx
j3 = w1[2] *w2[2] * wave_nm[2]^4/xx
.....
J_total =j1+j2+J3

print j_total
end

As you can see, I don't know how to introduce wave scaling into my parameters. I have 400 point to calculate ( from w1[0] to w1[400] )
W1 and w2 have already set wave scaling from 400 to 800.
wave_nm is another wave scale I have to introduce in order to take the equation in effect.
In other words, the display w1 vs wave_nm is the same as display w1

SO, i tried this
Function SpectralIntegral(w1,w2,wave_nm)
wave w1,w2,wave_nm
wavestats /q w2
variable xx = v_sum
variable i, j
for (i =0; i <10; i  +=1)
   
j = w1[i] *w2[i] * wave_nm[i]^4/xx
print j
endfor
end


I can get all the J but in print without integral. How can I add all the value from i<400?

This is my understanding of using IGOR. Can you give me suggestion for making the formula better?

August 24, 2015 at 05:12 pm - Permalink

In your second code example change

j = ....

to j += ...


This will accumulate the sum that you desire.


Function SpectralIntegral(w1,w2,wave_nm)
wave w1,w2,wave_nm
wavestats /q w2
variable xx = v_sum
variable i, j
for (i =0; i <400; i  +=1)
   
j += w1[i] *w2[i] * wave_nm[i]^4/xx

endfor
end



You could also take advantage of Igor's wave math capability and your wave scaling...

Function SpectralIntegral(w1,w2)
wave w1,w2
wavestats /q w2
variable xx = v_sum
duplicate w1, wProduct
variable vSum
wProduct() = w1(x) *w2(x) * (x)^4/xx
vSum = sum(wProduct)
print vSum
end


I haven't tested this code, so there may be and error here. Hopefully this will help, regardless.

August 25, 2015 at 05:39 am - Permalink

Thank you jtigor

The  += works

The  wProduct() = w1(x) *w2(x) * (x)^4/xx does not work.

A window with title "Function Compilation Error"

Error in
 wProduct() = w1(x) *w2(x) * (x)^4/xx

Saying:
Can't use wave(x)= in a function. Use x2pnt and wave[p]= instead

Any ideas how to alter the function of wave(x)?

Thanks!

August 25, 2015 at 09:22 am - Permalink

viralvector wrote:
Thank you for the input. I have to admit that I cannot figure out how to do it still.


It may be a good idea to take a couple of steps back to make sure that you understand what you are doing.

Suppose you have a wave representing the simple function f(x)=x. Your wave of 100 points starts at x=1 and goes to x=5.
make/N=100 wfx
SetScale/I x 1,5,"", wfx
wfx=x
display wfx


Now suppose you wanted to integrate this wave from x=2 to x=4. You can add the following two functions to your procedure window:
// this is the actual integrand.  It calls on a function that interpolates your wave data.
Function myUserFunction(inX)
    Variable inX
   
    Variable res=myFx(inX)
    return res
End

// this is a very simple function that interpolates wave data.  In a general case you can replace it with other types of interpolation.
Function myFx(inX)
    Variable inX
   
    Wave w=root:wfx
    return w(inx)
End


Clearly you can simplify the above but the structure I am suggesting is useful, e.g., if your interpolation involves XY pair of waves... Once you have these functions in place, you can simply execute:
print Integrate1D(myUserFunction,2,4,1)


This example should illustrate to you pretty much everything you need in order to obtain the integrals in question. I recommend that you do not attempt to sum up your data instead of using Integrate1D.

A.G.
WaveMetrics, Inc.

August 25, 2015 at 10:12 am - Permalink

viralvector wrote:


The  wProduct() = w1(x) *w2(x) * (x)^4/xx does not work.

A window with title "Function Compilation Error"

Error in
 wProduct() = w1(x) *w2(x) * (x)^4/xx

Saying:
Can't use wave(x)= in a function. Use x2pnt and wave[p]= instead

Any ideas how to alter the function of wave(x)?

Thanks!


Sorry, I was wrong. You cannot use the wave scaling like this in a function. You will need to use point scaling:

Function SpectralIntegral(w1,w2,wave_nm)
wave w1,w2,wave_nm
wavestats /q w2
variable xx = v_sum
duplicate w1, wProduct
variable vSum
wProduct[] = w1[p] * w2[p] * wave_nm[p]^4/xx
vSum = sum(wProduct)
print vSum
end



August 25, 2015 at 12:28 pm - Permalink