[Python-ideas] real numbers with SI scale factors

[Ken Kundert]

Chris,
    I was not able to get an astrophyics example, but I do have a reasonable one 
that performs a spectral analysis of the output of an analog to digital 
converter, something radio astronomers are known to do.

I am including the code, but it requires a rather large data file to run, which 
I will not include.  The code uses my 'engfmt' library from pypi to perform 
conversion to SI form. In this example, there is no need for conversion from SI 
form.

    #!/usr/bin/env python3

    import numpy as np
    from numpy.fft import fft, fftfreq, fftshift
    import matplotlib as mpl
    mpl.use('SVG')
    from matplotlib.ticker import FuncFormatter
    import matplotlib.pyplot as pl
    from engfmt import Quantity, set_preferences
    set_preferences(spacer=' ')

    def mag(spectrum):
        return np.absolute(spectrum)

    def freq_fmt(val, pos):
        return Quantity(val, 'Hz').to_eng()

    def volt_fmt(val, pos):
        return Quantity(val, 'V').to_eng()

    freq_formatter = FuncFormatter(freq_fmt)
    volt_formatter = FuncFormatter(volt_fmt)

    data = np.fromfile('delta-sigma.smpl', sep=' ')

    time, wave = data.reshape((2, len(data)//2), order='F')
    timestep = time[1] - time[0]
    nonperiodicity = wave[-1] - wave[0]
    period = timestep * len(time)
    print('timestep = {}'.format(Quantity(timestep, 's')))
    print('nonperiodicity = {}'.format(Quantity(nonperiodicity, 'V')))
    print('timepoints = {}'.format(len(time)))
    print('freq resolution = {}'.format(Quantity(1/period, 'Hz')))

    window = np.kaiser(len(time), 11)/0.37
        # beta=11 corresponds to alpha=3.5 (beta = pi*alpha), /.4 # 0.4 is the 
        # processing gain with alpha=3.5 is 0.37
    #window = 1
    windowed = window*wave

    spectrum = 2*fftshift(fft(windowed))/len(time)
    freq = fftshift(fftfreq(len(wave), timestep))

    fig = pl.figure()
    ax = fig.add_subplot(111)
    ax.plot(freq, mag(spectrum))
    ax.set_yscale('log')
    ax.xaxis.set_major_formatter(freq_formatter)
    ax.yaxis.set_major_formatter(volt_formatter)
    pl.savefig('spectrum.svg')
    ax.set_xlim((0, 1e6))
    pl.savefig('spectrum-zoomed.svg')

When run, this program prints the following diagnostics to stdout:

    timestep = 20 ns
    nonperiodicity = 2.3 pV
    timepoints = 27994
    freq resolution = 1.7861 kHz

It also generates two SVG files. I have converted one to PNG and attached it.

A few comments:
1. The data in the input file ('delta-sigma.smpl') has low dynamic range and is 
   machine generated, and not really meant for direct human consumption. As 
   such, it does not benefit from using SI scale factors.  But there are 
   certainly cases where the data has both high dynamic range and is intended 
   for people to examine it directly. In those cases it would be very helpful if 
   NumPy was able to directly read the file.  As the language exists today, 
   I would need to read the file myself, manually convert it, and feed the 
   result to NumPy.
2. Many of these numbers that are output do have high dynamic range and are 
   intended to be consumed directly by humans. These benefit from using SI scale 
   factors.  For example, the 'freq resolution' can vary from Hz to MHz and 
   'nonperiodicity' can vary from fV to mV.
3. Extra effort was expended to make the axis labels on the graph use SI scale 
   factors so as to make the results 'publication quality'.

My hope is that if Python accepted SI literals directly, then both NumPy nad 
MatPlotLib would also be extended to accept/use these formats directly, 
eliminating the need for me to do the conversions and manage the axes.

-Ken



On Mon, Aug 29, 2016 at 06:02:29PM +1000, Chris Angelico wrote:
> On Mon, Aug 29, 2016 at 5:07 PM, Ken Kundert
> <python-ideas at shalmirane.com> wrote:
> >
> > I talked to astrophysicist about your comments, and what she said was:
> > 1. She would love it if Python had built in support for real numbers with SI
> >    scale factors
> > 2. I told her about my library for reading and writing numbers with SI scale
> >    factors, and she was much less enthusiastic because using it would require
> >    convincing the rest of the group, which would be too much effort.
> > 3. She was amused by the "kilo pico speed of light" comment, but she was adamant
> >    that the fact that you, or some system administrator, does not understand
> >    what kpc means has absolutely no affect on her desired to use SI scale
> >    factors. Her comment: I did not write it for him.
> > 4. She pointed out that the software she writes and uses is intended either for
> >    herself of other astrophysicists. No system administrators involved.
> 
> So can you share some of her code, and show how the ability to scale
> unitless numbers would help it?
> 
> ChrisA
> _______________________________________________
> Python-ideas mailing list
> Python-ideas at python.org
> https://mail.python.org/mailman/listinfo/python-ideas
> Code of Conduct: http://python.org/psf/codeofconduct/
-------------- next part --------------
A non-text attachment was scrubbed...
Name: spectrum-zoomed.svg.png
Type: image/png
Size: 49011 bytes
Desc: not available
URL: <http://mail.python.org/pipermail/python-ideas/attachments/20160829/7dfd31cb/attachment-0001.png>