A simple model for measuring the period of a set of flame bursts: http://i.imgur.com/tbOM4dG.png

check_fft.py

# check_fft -- a simple little test of fft
import pylab
import scipy.special
import numpy as np
from pysurvey import plot as pplot



def _gauss(x, a, x0, sigma):
    '''Commented out the gaussian and now a skewed gaussian.'''
    # return a*np.exp((-(x-x0)**2/float(sigma)))
    alpha = 4.0 
    tmp = 0.5*(1+scipy.special.erf(alpha*(x-x0)))
    return a*np.exp((-(x-x0)**2/float(sigma)))*tmp


def main():
    
    # Nominal frequencies for the two signals
    freq1 = 1/30.
    freq2 = 1/20.
    
    # take measurements every 0.1s for 1000s
    delta = 0.1
    t = np.arange(0,1000,delta)
    v = np.zeros(len(t))
    
    # add a normal spike every 10 sec with amplitude == 1
    # the normal spikes have width 2.0
    # Do not cover the entire range so that you can see that it is not
    # just periodicity that is changing things.
    tmp = np.arange(10,400,1/freq1)
    for tt in tmp:
        v += _gauss(t, 2.0, tt, 5.0)
    
    # have a big spike every 50s with amplitude == 5
    # since they are larger have a larger width of 4
    # start the big ones at 30
    # add some jitter(sigma=2.0) to the start times
    tmp = np.arange(30,200,1/freq2) 
    tmp2 = np.random.normal(0,2.0, len(tmp))
    for tt,jit in zip(tmp,tmp2):
        v += _gauss(t, 5.0, tt+jit, 5.0)
    
    # add some noise -- because this is the real world
    # signal to noise of 0.2 for the small ones
    v += np.random.normal(0,0.2, len(v))
    
    
    pylab.figure(1, figsize=(12,12))
    
    # Show the peak model
    pplot.setup(subplt=(2,2,1), xlabel='Time[s]', ylabel='Amplitude[AU:Arbitrary Units]', 
                subtitle='Skew Normal Model')
    x = np.arange(0,10,0.1)
    pylab.plot(x,_gauss(x,1, 5, 2))
    # pplot.line(x=5, color='orange')
    
    
    # show the time series
    pplot.setup(subplt=(2,2,2), 
                xlabel='Time[s]', 
                ylabel='Amplitude[AU:Arbitrary Units]',
                subtitle='Time Series')
    pylab.plot(t,v)
    
    # Calculate fft
    freq = np.fft.fftfreq(t.size, d=delta)
    sp = np.fft.fft(v)
    ps = np.abs(sp)**2.0

    # Show the power spectrum
    pplot.setup(subplt=(2,1,2), 
                xr=[0,0.2], xlabel='Frequency[Hz]',
                ylabel='Power Spectrum',
                subtitle='Nominal frequencies at {:0.2f} and {:0.2f}\nFirst order harmonics shown'.format(freq1,freq2))
    idx = np.argsort(freq)
    pylab.plot(freq[idx], ps[idx])
    pplot.line(x=[freq1,freq2], color='orange')
    pplot.line(x=[2*freq1,2*freq2], color='orange', linestyle='--')
    pylab.show()
    
    
    
    
    
    
    








if __name__ == '__main__':
    main()