On Tue, Mar 31, 2009 at 8:54 PM, Jochen S <[email protected]> wrote:
> On Tue, Mar 31, 2009 at 7:13 AM, João Luís Silva <[email protected]> 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 >> [email protected] >> http://mail.scipy.org/mailman/listinfo/numpy-discussion >> > >
_______________________________________________ Numpy-discussion mailing list [email protected] http://mail.scipy.org/mailman/listinfo/numpy-discussion
