Welcome, guest | Sign In | My Account | Store | Cart

Real_FFT wraps the gsl_fft_real_transform in a python3 setting. This recipe serves as a complete example for "Recipe 576549: gsl with python3".

Python, 104 lines
  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
'''
    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]

gsl_setup.py is the code of the "gsl with python3" recipe.

A Real_FFT object is callable. Called with a real array it returns the FFT. The result is re-arranged from the gsl order to python complex type from low to high frequency. The result does not include the negative frequencies since the amplitudes are symmetric. Other methods in the class are produce_frequency and produce_period to give frequency and period in your units. The samples keyword argument describes the length of the transformed array. It defaults to the length of the previous FFT computed by the object. produce_frequency and produce_period accept at most one of the keyword only arguments sample_interval, sample_rate, or total_length. The default is sample_interval = 1 project dependent unit.

>>> ystar = real_fft(y)
>>> T = real_fft.produce_period()
>>> f = real_fft.produce_frequency()
>>> # now ystar[i] corresponds to f[i] and T[i]

The object caches the previous precomputed wavetable for reuse with same length transform. Cache size == 1. Real_FFT supports del so that memory can be easily reclaimed. This behavior differs from that of scipy which can (and on my computer did) run out of memory for transforms on many different lengths.

1 comment

David Lambert (author) 13 years, 1 month ago  # | flag

I've wrapped the high level gsl_interp routines and the outstanding gsl_rng routines for python3k. By request, I'll post these also.