gsl real fft in python3 (Python recipe) by David Lambert

ActiveState Code (http://code.activestate.com/recipes/576550/)

'''
    provide a simple python3 interface to the gsl_fft_real_transform function
'''

import sys
import itertools
from gsl_setup import *

def grouper(n, iterable, fillvalue=None):
    # http://docs.python.org/dev/3.0/library/itertools.html#module-itertools
    "grouper(3, 'ABCDEFG', 'x') --> ABC DEF Gxx"
    args = [iter(iterable)] * n
    return itertools.zip_longest(fillvalue=fillvalue, *args)

real_workspace_alloc = setup(
    gsl.gsl_fft_real_workspace_alloc,[c_ulong,],c_void_p)
real_wavetable_alloc = setup(
    gsl.gsl_fft_real_wavetable_alloc,[c_ulong,],c_void_p)
real_workspace_free =setup(gsl.gsl_fft_real_workspace_free ,[c_void_p,])
real_wavetable_free =setup(gsl.gsl_fft_real_wavetable_free ,[c_void_p,])

real_transform = setup(gsl.gsl_fft_real_transform,
                       [c_void_p,c_ulong,c_ulong,c_void_p,c_void_p],)


class Real_FFT:

    '''
        returns the complex values of the real transform of the real data.
        return value[0] describes the offset,
                    [1] is amplitude of term for wavelength = data length
                    etceteras
                    [-1] amp of wavelength = twice sample distance
    '''

    def __init__(self):
        self.n = 0

    def __call__(self,data):
        if len(data) < 2:
            if 1 == len(data):
                return data[:]
            return []
        if len(data) != self.n:
            self.__del__()
            self.n = len(data)
            size = c_ulong(self.n)
            self.workspace = real_workspace_alloc(size)
            self.wavetable = real_wavetable_alloc(size)
        a = array('d',data)       # need a copy of the data
        real_transform(ADDRESS(a),1,self.n,self.wavetable,self.workspace)
        rv = [complex(a[0]),]
        rv.extend(itertools.starmap(complex,grouper(2,a[1:],fillvalue=0)))
        return rv

    def __del__(self):
        if self.n:
            try:
                real_workspace_free(self.workspace)
                real_wavetable_free(self.wavetable)
            except AttributeError:
                print('Attribute error while freeing FFT auxiliary storage',
                      file=sys.stderr)
            except:
                print('error freeing FFT auxiliary storage',
                      file=sys.stderr)

    def produce_frequency(self,*,samples=None,sample_interval=None,sample_rate=None,total_length=None):
        '''
            return the frequency grid based on actual sizes (default sample_interval=1).
        '''
        n = samples or self.n
        if not n:
            return array('d')
        args_specified = 3 - ((not sample_interval)+(not sample_rate)+(not total_length))
        if 1 < args_specified:
            raise TypeError('specify at most one of [sample_rate, total_length, sample_interval]')
        if 0 == args_specified:
            L = n
        elif sample_interval:
            L = n*sample_interval
        elif sample_rate:
            L = n/sample_rate
        else:
            L = total_length
        return as_array(waves/L for waves in range(1+n//2))
            
    def produce_period(self,*args,**kwargs):
        '''
            return the period grid based on actual sizes.
            frequency of zero --> period 0.  what else to do?
        '''
        f2T = self.produce_frequency(*args,**kwargs)
        for i in range(1,len(f2T)):
            f2T[i] = 1/f2T[i]
        return f2T

real_fft = Real_FFT()

def magnitude(a):
    return [abs(b) for b in a]

def phase(a):
    return [phase(b) for b in a]