[Numpy-discussion] Recognizing a cycle in a vector

Manolo Martínez manolo at austrohungaro.com
Thu Dec 3 03:45:57 EST 2015


> > >> Is there any way to check for cycles in this situation?
> > 
> > > Fast fourier transform (fft)?
> > 
> > +1 For using a discrete Fourier transform, as implemented by numpy.fft.fft.
> > You mentioned that you sample at points which do not correspond with the
> > period of the signal; this introduces a slight complexity in how the
> > Fourier transform reflects information about the original signal. I attach
> > two documents to this email with details about those (and other)
> > complexities. There is also much information on this topic online and in
> > signal processing books.
> 

So, I thought I'd report back on what I've ended up doing. Given that
the cycles I find in my data are usually very close to sine waves, the
following works well enough:


    def periodic_vector(vector):
        """
        Take the FFT of a vector, and eliminate all components but the
        two main ones (i.e., the static and biggest sine amplitude) and
        compare the reconstructed wave with the original. Return true if
        close enough 
        """
        rfft = np.fft.rfft(vector)
        magnitudes = np.abs(np.real(rfft))
        choice = magnitudes > sorted(magnitudes)[-3]
        newrfft = np.choose(choice, (np.zeros_like(rfft), rfft))
        newvector = np.fft.irfft(newrfft)
        return np.allclose(vector, newvector, atol=1e-2)


This is doing the job for me at the moment, but there are, that I can
see, a couple of things that could be improved (and surely more that I
cannot see): 

1) this func sorts the absolute value of the amplitudes to find the two
most important  components, and this seems overkill for large vectors. 

2) I'm running the inverse FFT, and comparing to the initial vector, but
it should be possible to make a decision solely based on the size of
terms in the FFT itself. I'm just not confident enough to design a test
based on that.

Anyway, thanks to those who pointed me in the right direction.

Manolo



More information about the NumPy-Discussion mailing list