[Numpy-discussion] Optical autocorrelation calculated with numpy is slow

[Jochen S]

On Tue, Mar 31, 2009 at 8:54 PM, Jochen S <cycomanic at gmail.com> wrote:

> On Tue, Mar 31, 2009 at 7:13 AM, João Luís Silva <jsilva at fc.up.pt> wrote:
>
>> Hi,
>>
>
>
>> I wrote a script to calculate the *optical* autocorrelation of an
>> electric field. It's like the autocorrelation, but sums the fields
>> instead of multiplying them. I'm calculating
>>
>> I(tau) = integral( abs(E(t)+E(t-tau))**2,t=-inf..inf)
>>
>
> An autocorrelation is just a convolution, which is a multiplication in
> frequency space. Thus you can do:
> FT_E = fft(E)
> FT_ac=FT_E*FT_E.conj()
> ac = fftshift(ifft(FT_ac))
>
> where E is your field and ac is your autocorrelation. Also what sort of
> autocorrelation are you talking about. For instance SHG autocorrelation is
> an intensity autocorrelation thus the first line should be:
> FT_E = fft(abs(E)**2)
>

Sorry I was reading over your example to quickly earlier, you're obviously
using intensity autocorrelation so what you should be doing is:
FT_E=fft(abs(E)**2)
FT_ac = FT_E*FT_E.conj()
ac = fftshift(ifft(FT_ac))



>
> HTH
> Jochen
>
>
>> with script appended at the end. It's too slow for my purposes (takes ~5
>> seconds, and scales ~O(N**2)). numpy's correlate is fast enough, but
>> isn't what I need as it multiplies instead of add the fields. Could you
>> help me get this script to run faster (without having to write it in
>> another programming language) ?
>>
>> Thanks,
>> João Silva
>>
>> #--------------------------------------------------------
>>
>> import numpy as np
>> #import matplotlib.pyplot as plt
>>
>> n = 2**12
>> n_autocorr = 3*n-2
>>
>> c = 3E2
>> w0 = 2.0*np.pi*c/800.0
>> t_max = 100.0
>> t = np.linspace(-t_max/2.0,t_max/2.0,n)
>>
>> E = np.exp(-(t/10.0)**2)*np.exp(1j*w0*t)    #Electric field
>>
>> dt = t[1]-t[0]
>> t_autocorr=np.linspace(-dt*n_autocorr/2.0,dt*n_autocorr/2.0,n_autocorr)
>> E1 = np.zeros(n_autocorr,dtype=E.dtype)
>> E2 = np.zeros(n_autocorr,dtype=E.dtype)
>> Ac = np.zeros(n_autocorr,dtype=np.float64)
>>
>> E2[n-1:n-1+n] = E[:]
>>
>> for i in range(2*n-2):
>>     E1[:] = 0.0
>>     E1[i:i+n] = E[:]
>>
>>     Ac[i] = np.sum(np.abs(E1+E2)**2)
>>
>> Ac *= dt
>>
>> #plt.plot(t_autocorr,Ac)
>> #plt.show()
>>
>> #--------------------------------------------------------
>>
>> _______________________________________________
>> Numpy-discussion mailing list
>> Numpy-discussion at scipy.org
>> http://mail.scipy.org/mailman/listinfo/numpy-discussion
>>
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.python.org/pipermail/numpy-discussion/attachments/20090331/d5771287/attachment.html>