Charles R Harris wrote:

That should work. The first two integrals are actually the same, but need to be E(t)*E(t).conj(). The second integral needs twice the real part of E(t)*E(t-tau).conj(). Numpy correlate should really have the conjugate built in, but it doesn't.

Chuck

It worked, thanks. João Silva

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

On Tue, Mar 31, 2009 at 7:13 AM, João Luís Silva <jsilva@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@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion