Re: [Numpy-discussion] numpy.correlate with phase offset 1D data series
Thank you for the input! It sounds like Fourier methods will be fastest, by design, for sample counts of hundreds to thousands. I currently do steps like: Im1 = get_stream_array_data() Im2 = load_template_array_data(fh2) ##note: len(im1)==len(im2) Ffft_im1=fftpack.rfft(Im1) Ffft_im2=fftpack.rfft(Im2) R1= (Ffft_im1 * Ffft_im2.conjugate()) R2= (abs(Ffft_im1) * abs(Ffft_im2)) R = R1 / R2 IR=fftpack.irfft(R) flat_IR = numpy.ravel(numpy.transpose(IR)).real I= numpy.argmax(flat_IR) phase_offset = (I % len(Im1)) At 09:29 AM 3/4/2008, Anne Archibald wrote: * What do you want to happen at the endpoints? Without padding, only a small interval (the difference in lengths plus one) is valid. Zero-padding works, but guarantees a fall-off at the ends. Circular correlation is easy to implement but not appropriate most of the time. How much should I be concerned?, since the only desired information from this is the scalar best-fit phase value, presumably the argmax() of the xcorr. In current operation, imagine a tone pattern/template of n samples which we want to align to streaming data; the desired result (at least in my current FFT code) is the sample number of recent ADC data where the zero'th sample of the pattern best aligns. Since it is a repeating pattern, we know that it will always align somewhere in the latest n samples. * Do you care about sub-sample alignment? How much accuracy do you really need? Integer alignment is sufficient, due both to electronic noise, and desired phase The other common application is to have a template (that presumably falls to zero at its endpoint) and to want to compute a running correlation against a stream of data. This too can be done both ways, depending on the size of the template; all that is needed is to think carefully about overlaps. This is very much what the application is, although the template does not terminate at zero. It does terminate at a value near the zero'th value however, and I assumed the FFTs would be well-behaved. Ray -- No virus found in this outgoing message. Checked by AVG Free Edition. Version: 7.5.516 / Virus Database: 269.21.4/1310 - Release Date: 3/4/2008 8:35 AM ___ Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] numpy.correlate with phase offset 1D data series
At 03:28 PM 3/3/2008, Ann wrote: Sounds familiar. If you have a good signal-to-noise ratio, you can get subpixel accuracy by oversampling the irfft, or better but slower, by using numerical optimization to refine the peak you found with argmax. the S/N here is poor, and high data rates work against me too... I would be inclined to use convolve (or scipy.ndimage.convolve, which uses a Fourier-domain method), since it is somewhat better specified. I'll give it a try as well. I'm guessing scipy.ndimage.correlate1d is a Fourier method too? From: Timothy Hochberg [EMAIL PROTECTED] I'm going to guess that you are using some flavor of Unix, since I also downloaded using Firefox and the data ends up corrupted. My hypothesis is that Firefox doesn't recognize the mime type and treats it as a text file, corrupting it on Windows, but not on Unix. Then again, maybe you're not using Unix and my installation of Firefox is just broken. I think that is the case, I have Win2K on this box With the CSV version I do get a peak at the (un)expected location (7489//2). The peak is pretty flat and only twice the size of the surrounding gunk, but it looks more or less legit. I don't see that in my pylab plot! There's actually a dip, and the whole plot is symmetric about 3744 http://rjs.org/Python/corr_array.jpg, self xcorr of http://rjs.org/Python/data.jpg I'll be upgrading my install here shortly though to py2.5 and associated libs. My compiler/distutils environment is broken. Thanks, Ray -- No virus found in this outgoing message. Checked by AVG Free Edition. Version: 7.5.516 / Virus Database: 269.21.4/1310 - Release Date: 3/4/2008 8:35 AM ___ Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion
[Numpy-discussion] numpy.correlate with phase offset 1D data series
I'm trying to figure out what numpy.correlate does, and, what are people using to calculate the phase shift of 1D signals? (I coded on routine that uses rfft, conjugate, ratio, irfft, and argmax based on a paper by Hongjie Xie An IDL/ENVI implementation of the FFT Based Algorithm for Automatic Image Registration - but that seems more intensive than it could be.) In numpy, an identity import numpy from pylab import * l=[1,5,3,8,15,6,7,7,9,10,4] c=numpy.correlate(l,l, mode='same') plot(c) peaks at the center, x=5, and is symmetric when the data is rotated by 2 c=numpy.correlate(l, l[-2:]+l[:-2], mode='same') it peaks at x=3 I was expecting, I guess, that the peak should reflect the x axis shift, as in http://en.wikipedia.org/wiki/Cross-correlation#Explanation If I use a real time domain signal like http://rjs.org/Python/sample.sig fh = open(r'sample.sig','rb') s1 = numpy.fromstring(fh.read(), numpy.int32) fh.close() an identity like c=numpy.correlate(s1, s1, mode='same') plots like noise. No virus found in this outgoing message. Checked by AVG Free Edition. Version: 7.5.516 / Virus Database: 269.21.3/1308 - Release Date: 3/3/2008 10:01 AM ___ Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion
[Numpy-discussion] numpy.correlate with phase offset 1D data series
At 01:24 PM 3/3/2008, you wrote: If you use 'same' or 'full' you'll end of with different amounts of offset. I imagine that this is due to the way the data is padded. The offset should be deterministic based on the mode and the size of the data, so it should be straightforward to compensate for. I agree If I use a real time domain signal like http://rjs.org/Python/sample.sig fh = open(r'sample.sig','rb') s1 = numpy.fromstring(fh.read(), numpy.int32) fh.close() When I download this, it's full of NaNs. There's either a problem in the way I downloaded it or in the uploaded file. You didn't by chance upload it as an ASCII file did you? I just tested the URL myself with Firefox; it came down OK. It is a binary string from numpy.tostring(), 29,956 bytes of int32. It has a fundamental of 42 cycles in the data, and other fs of less power. I just uploaded a http://rjs.org/Python/sample.csv version Xie's 2D algorithm reduced to 1D works nicely for computing the relative phase, but is it the fastest way? It might be, since some correlation algorithms use FFTs as well. What does _correlateND use, in scipy? Thanks, Ray -- No virus found in this outgoing message. Checked by AVG Free Edition. Version: 7.5.516 / Virus Database: 269.21.3/1308 - Release Date: 3/3/2008 10:01 AM ___ Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion
[Numpy-discussion] image to array doubt
At 10:00 AM 2/29/2008, you wrote: Robin wrote I'm not sure why they would be doing this - to me it looks they might be using Image as a convenient way to store some other kind of data... thanks Robin, I am wondering if there is a more straightforward way to do these.. especially the vector to image function I would normally suggest PIL, as it contains most all of these types of arrayimage functions See: http://www.pythonware.com/library/pil/handbook/imagemath.htm http://www.pythonware.com/library/pil/handbook/imageops.htm etc... as well as the built-ins for numpy array objects http://effbot.org/zone/pil-changes-116.htm frombuffer, fromstring, fromarray, tostring etc. http://www.pythonware.com/library/pil/handbook/image.htm (I've used them for home astronomy projects, myself.) Ray Schumacher Blue Cove Interactive No virus found in this outgoing message. Checked by AVG Free Edition. Version: 7.5.516 / Virus Database: 269.21.2/1304 - Release Date: 2/29/2008 8:18 AM ___ Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion
[Numpy-discussion] numpy Intel ICC build help?
We've just built Python 2.4 in the lab with Intel ICC and MS VS 2005, but had problems building 2.5, and also numpy, with the MKL. Would someone be willing to share their build experience or project files? Ray -- No virus found in this outgoing message. Checked by AVG Free Edition. Version: 7.5.516 / Virus Database: 269.19.2/1222 - Release Date: 1/13/2008 12:23 PM ___ Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] parallel numpy - any info?
At 04:27 AM 1/8/2008, you wrote: 4. Re: parallel numpy (by Brian Granger) - any info? (Matthieu Brucher) From: Matthieu Brucher [EMAIL PROTECTED] MKL does the multithreading on its own for level 3 BLAS instructions (OpenMP). There was brief debate yesterday among the Pythonians in the lab as to whether the MKL operates outside of the GIL, but it was general understanding that it does. It is still unclear to me whether Python/numpy compiled with MKL would be freely re-distributable, as the MSVC version is. Then you must use tools like the processing module, MPI, ... I have been using shared memory and numpy with good results lately on Win32 with tagged memory. http://projects.scipy.org/pipermail/numpy-discussion/2007-October/029505.html On *NIX, arrayfrombuffer(), mmap and MAP_SHARED, or POSH http://poshmodule.sourceforge.net/ seem like the ticket. Whatever happened to pyGA? GA http://www.emsl.pnl.gov/docs/global/ is still around. http://www.ece.lsu.edu/jxr/pohll-02/papers/jarek.pdf Best, Ray Schumacher -- No virus found in this outgoing message. Checked by AVG Free Edition. Version: 7.5.516 / Virus Database: 269.17.13/1211 - Release Date: 1/6/2008 11:57 AM ___ Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Changing the distributed binary for numpy 1.0.4 for windows ?
At 02:32 AM 12/11/2007, you wrote: If so I'd be happy to contribute part of the purchase price, and I assume others would too. What's more, I *have* an old PIII at home. The main company I consult for is set to buy the Intel compiler and FFT lib for Windows, for the express purpose of compiling Python, numpy, and the fastest FFT for each CPU, for new products. I develop on both an old dual PIII and a new Core Duo, and there is also every other flavor in the shop. As I read it, the binaries should be distributable. The catch is, I have never compiled Python or numpy - I had seen a poster who offered help. When the company actually purchases the product I'd be glad to do it on 2-3 targets if someone can assist with the parameters. We have one consultant here who has done it on Linux. Ray Schumacher Congitive Vision 8580 Production Ave., Suite B San Diego, CA 92121 858.578.2778 http://cognitivevision.com/ -- No virus found in this outgoing message. Checked by AVG Free Edition. Version: 7.5.503 / Virus Database: 269.17.0/1180 - Release Date: 12/10/2007 2:51 PM ___ Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] numpy FFT memory accumulation
At 10:57 PM 11/1/2007, Charles R Harris wrote: An additional complication is that I pass the numpy (or Numeric) array address to the ctypes library call so that the data is placed directly into the array from the call. I use the if/else end wrap logic to determine whether I need to do a split and copy if the new data wraps. OK. Hmm, I wonder if you would lose much by taking a straight forward radix-2 fft and teaching it to use modular indices? Probably not worth the trouble, but an fft tailored to a ring buffer might be useful for other things. The problem is, I once compiled my own FFT dll to call from Python and it was 2x slower than FFTPACK - I'm not math-smart enough to make all of the caching and numerical shortcuts. That, and Intel's optimized FFT is 3x faster than FFTPACK... I may still try to do a zoom/range FFT which does not compute all bins, and maybe only with a sine transform, which (I think) should be sufficient to determine peak real frequency and maybe use fewer cycles. Probably the easiest thing is to just copy the ring buffer out into a linear array. I do that for the buffer-wrap condition, while simply assigning a slice (no copy) to the temp array otherwise. I used Numeric functions for the ~40% speed increase, but I don't I know that numarray was slow in creating small arrays, but is Numpy really that bad compared to Numeric? I just saw the # of FFTs/sec go from 390 to 550 just by switching numpy to Numeric (Intel Core Duo). Add a timer to my previous code posts and see how your results look. For the mega-arrays a lot of the numpy developers work with it is much faster, and I now find Numeric is lacking many other niceties, like frombuffer(). Ray -- No virus found in this outgoing message. Checked by AVG Free Edition. Version: 7.5.503 / Virus Database: 269.15.18/1104 - Release Date: 11/1/2007 6:47 PM ___ Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] numpy FFT memory accumulation
At 11:55 PM 10/31/2007, Travis wrote: Ray S wrote: I am using fftRes = abs(fft.rfft(data_array[end-2**15:end])) At first glance, I would say that I don't expect memory to be growing here, so it looks like a problem with rfft that deserves looking into. I saw that Numeric did also (I still use Numeric for smaller array speed) but much more slowly. I will try to repeat with a small demo and post. Does data_array keep growing? no, it is a 64k circular buffer Which reminds me, I've been wanting to try to build a Numeric circular buffer object; one that, when sliced or assigned to, auto-magically wraps as needed, exposing only an extra pointer attribute for the end. I currently always do it with Python if-s and pointer vars for these products that only do real time data analysis. From: Anne Archibald [EMAIL PROTECTED] If the range is *really* small, you can try using a DFT - sometimes that is fast enough, and gives you just the bins you're curious about. I've considered weave'ing a simple sine transform with specified range, but until I do and test I won't know if my own implementation is any faster than just the FFTPACK. If the range is bigger than that, but still a small fraction of the FFT size, you can do some tricks where you band-pass filter the data ... There are also zoom fft and chirp-z techniques which are supposed to give you only part of the FFT, but the wisdom is that unless you want less than a few percent of the data points you're better just FFTing and throwing lots of data away. I've tried zoom before; the issue was just that - 2 FFTs and a shift or convolution eats a lot of CPU cycles and falls behind the real time data. The range of interest in the Fourrier domain is small, 3kHz-7kHz. The sample rate is high for precise phase information. I've got some more testing to do, it seems. Thanks, Ray -- No virus found in this outgoing message. Checked by AVG Free Edition. Version: 7.5.503 / Virus Database: 269.15.15/1101 - Release Date: 10/31/2007 10:06 AM ___ Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] numpy/Windows shared arrays between processes?
At 05:22 AM 10/9/2007, David Cournapeau wrote: Could not this be because you compiled the posh sources with a compiler/runtime which is different than the other extensions and python interpreter ? It definitely was - since my 2.4 wanted the free 7.1 compiler, I (and anyone else who didn't download it in time) are now seemingly SOL since it is no longer available. I saw much discussion of this as well, but even 2.5 is now fixed on 7.1 and reports of compiling distutil modules with the new MS SDK and having them work at all with 2.4 were very mixed. I also tried GCC and had a litany of other errors with the posh. Sebastian Haase added: I was in fact experimenting with this. The solution seemed to lie in simple memmap as it is implemented in Windows: snip I had just found and started to write some tests with that MS function. If I can truly write to the array in one process and instantly read it in the other I'll be happy. Did you find that locks or semaphores were needed? (( I have to mention, that I could crash a process while testing this ... )) That was one of my first results! I also found that using ctypes to create arrays from the other process's address and laying a numpy array on top was prone to that in experimentation. But I had the same issue as Mark Heslep http://aspn.activestate.com/ASPN/Mail/Message/ctypes-users/3192422 of creating a numpy array from a raw address (not a c_array). Thanks, Ray Schumacher ___ Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Compile extension modules with Visual Studio 2005
Geoffrey Zhu wrote: Hi, I am about to write a C extension module. C functions in the module will take and return numpy arrays. I found a tutorial online, but I am not sure about the following: I agree with others that ctypes might be your best path. The codeGenerator is magic, if you ask me: http://starship.python.net/crew/theller/ctypes/old/codegen.html But, if the function is simple, why not weave.inline? What I have done is run the function once, hunt down the long-named library, copy it to the local directory, then include it explicitly and call its function. This eliminates some overhead time for the call. I use it to convert packed IEEE data from an ADC data read function, and it's faster than the manufacturer's own function version that returns scaled integers! Ray ___ Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] numpy array sharing between processes?
Andrew added: I'll pitch in a few donuts (and my eternal gratitude) for an example of shared memory use using numpy arrays that is cross platform, or at least works in linux, mac, and windows. I thought that getting the address from the buffer() of the array and creating a new one from it in the other process is most likely to work cross-platform, but I don't know enough of Python innards to speak confidently. I have used it with numarray before : N = numarray.zeros((1000,), numarray.Float) N.info() ... data: memory at 009d67b8 with size:8000 held by object 009d6798... I'm also considering trying more mix with ctypes, as in: http://www.scipy.org/Cookbook/Ctypes#head-db3c64e7ade3f0257084b9c03f03d286ff3b1b4f - Ray ___ Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] NumPy benchmark
I'm still curious about the licensing aspects of using Intel's compiler and libs. Is the compiled Python/numpy result distributable, like any other compiled program? Ray ___ Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] building numpy with Intel MS compiler and FFT
Thanks Rex, I'll give it a try next week. I've compiled both Numpy and Python 2.5 with the Intel compiler. On a Core 2 Duo, at least, the speed increase on Pybench was ~49%, even before compiling Python with icc. My post about it was on 25 Jan, and has subject: Compiling Python with icc ___ Numpy-discussion mailing list [EMAIL PROTECTED] http://projects.scipy.org/mailman/listinfo/numpy-discussion
