telegraphic / pfb_introduction

An interactive introduction to the polyphase filterbank technique for radio astronomy spectrometers

Polyphase filterbanks: an interactive introduction

This repository provides some useful materials for understanding polyphase filterbanks (PFBs). Foremostly, this repo provides an interactive notebook that gives a gentle walkthrough of how polyphase filterbanks work. Also provided is Python PFB code, some supplementary diagrams and a list of useful references.

Getting started

While you can view the notebook on github, it's more illustrative to run it youself. To run the notebook on your local computer, you'll need to clone this repo and have a working scientific python environment installed (jupyter, numpy, scipy and seaborn). If you don't have this set up, the Anaconda python environment is a good choice.

To run the notebook, start the jupyter notebook server running from the pfb_introduction directory, by running on the command line:

> jupyter notebook

This will open your browser, from which you can launch the pfb_introduction.ipynb notebook.

License

The materials herein are licensed under Creative Commons CC-BY 4.0.

Citing

If you find this useful, please consider citing the accompanying Chapter of the Handbook of Astronomical Instrumentation. This is available on arxiv:

Danny C. Price, Spectrometers and Polyphase Filterbanks in Radio Astronomy, 2016. arXiv 1607.03579

URL: http://arxiv.org/abs/1607.03579

pfb_introduction's Issues

Test using continuous uniform pdf instead of Gaussian distribution pdf

First of all good job with the code. I ported to Matlab and I could not find any error or problem.

The only issue that I have seen is when you test the PFB you use np.random.random() to generate noise. np.random.random has a uniform probability density function (pdf), according to the Numpy webpage:

Results are from the “continuous uniform” distribution over the stated interval. [0,1] by defeault

Thermal noise, which is what you want to simulate to emulate RF noise, has a Gaussian pdf. To simulate Gaussian noise you should use np.random.randn() instead, according to the Numpy webpage:

randn generates an array of shape (d0, d1, ..., dn), filled with random floats sampled from a univariate “normal” (Gaussian) distribution of mean 0 and variance 1

This code issue is in both the .py and the jupyter notebook.

Best Wishes!

Divide by zero in Testing leakage with sin waves [16]

At the end of running the notebook, I saw this:

<ipython-input-2-5aa2c85ecb9b>:12: RuntimeWarning: divide by zero encountered in log10
  return 10*np.log10(x)
<ipython-input-2-5aa2c85ecb9b>:12: RuntimeWarning: divide by zero encountered in log10
  return 10*np.log10(x)
<ipython-input-2-5aa2c85ecb9b>:12: RuntimeWarning: divide by zero encountered in log10
  return 10*np.log10(x)

Modified db() to find out why:

def db(x):
    """ Convert linear value to dB value """
    with np.errstate(divide='raise'):
        try:
            result = 10.0 * np.log10(x)
        except FloatingPointError as err:
            print("*** db: Oops, RunTimeWarning err:{}\nx={}".format(err, x))
            sys.exit(86)
    return result

*** db: Oops, RunTimeWarning err:divide by zero encountered in log10
x=[0.0, 0.44679333752998857, 0.44306315393081414, 0.4442646162251902, 0.440156033570239, 0.4413821996488543, 0.43824838511168346, 0.44033032545309336, 0.4385075516390453, 0.4407382696622541, 0.438755195229924, 0.43862168999098533, 0.43353566476580385, 0.4276080559521242, 0.4163105772054383, 0.4017930589462105, 0.3824773234985193, 0.358450440997461, 0.3313642399561077, 0.299138611002753, 0.2665950939027449, 0.2294763558673873, 0.19570095063646759, 0.15809497903660172, 0.13049071225171915, 0.10062253507273337, 0.0716498912407578, 0.0517037418281012, 0.03500554602925388, 0.02287041244613334, 0.01409740301219314, 0.008215325616891212, 0.004469257067938133, 0.0022346449331279443, 0.0010102222431275796, 0.0003952373260176991, 0.00012630368122208, 2.7618101912515074e-05, 2.127615890847873e-06, 4.860906435258243e-07, 2.8330137937249312e-06, 3.615253882009608e-06, 2.8720644548633765e-06, 1.673762482388919e-06, 7.627812308590531e-07, 2.7016636896519775e-07, 7.591345867972363e-08, 1.8882039831131855e-08, 6.461347548593791e-09, 5.169225787640778e-09, 4.7291286173112145e-09, 2.8045225869571647e-09, 9.93368827776587e-11, 3.738407450056794e-09, 2.242347297644185e-08, 5.63618320041525e-08, 9.454287791446423e-08, 1.1630420679175998e-07, 1.0932152533366832e-07, 7.390402146886183e-08, 2.9868511236445555e-08, 2.152625203499346e-09, 7.789024721578119e-09, 4.447370768285094e-08, 9.168721705976398e-08, 1.2327215881864552e-07, 1.2103213874988774e-07, 8.776368310067445e-08, 4.090986072926849e-08, 6.792643238846288e-09, 2.19668864711491e-09, 2.782163633755817e-08, 6.635125777562617e-08, 9.881463894052946e-08, 1.0351578526098045e-07, 7.680084200218828e-08, 4.002757625766943e-08, 9.199233876330439e-09, 3.19115890490225e-10, 1.5611848467412966e-08, 4.410871342345218e-08, 6.891131806253674e-08, 7.557326604980442e-08, 6.169749743287123e-08, 3.473783407151327e-08, 1.0090183910022984e-08, 7.828829157793151e-12, 8.407947966440945e-09, 2.8637561802230894e-08, 4.845193112634732e-08, 5.6149958709493355e-08, 4.833789301084396e-08, 2.9418708026228782e-08, 1.0113557978984861e-08, 3.087655437950417e-10, 4.241781579308114e-09, 1.8526256856559443e-08, 3.3786335472489205e-08, 4.519306841951586e-08, 3.881448924698381e-08, 2.4719796087206734e-08]

As you can see, the first element is 0 and it comes from here:

plt.plot(period, db(chan0_val)) # [16]

So, all 3 db() calculations fail with arguments chan0_val, chan1_val, and chan2_val. All 3 channel arrays have a 0 in the first element.

DEBUG chan0_val=[0.0, 0.44679333752998857, 0.44306315393081414, 0.4442646162251902, 0.440156033570239, 0.4413821996488543, 0.43824838511168346, 0.44033032545309336, 0.4385075516390453, 0.4407382696622541, 0.438755195229924, 0.43862168999098533, 0.43353566476580385, 0.4276080559521242, 0.4163105772054383, 0.4017930589462105, 0.3824773234985193, 0.358450440997461, 0.3313642399561077, 0.299138611002753, 0.2665950939027449, 0.2294763558673873, 0.19570095063646759, 0.15809497903660172, 0.13049071225171915, 0.10062253507273337, 0.0716498912407578, 0.0517037418281012, 0.03500554602925388, 0.02287041244613334, 0.01409740301219314, 0.008215325616891212, 0.004469257067938133, 0.0022346449331279443, 0.0010102222431275796, 0.0003952373260176991, 0.00012630368122208, 2.7618101912515074e-05, 2.127615890847873e-06, 4.860906435258243e-07, 2.8330137937249312e-06, 3.615253882009608e-06, 2.8720644548633765e-06, 1.673762482388919e-06, 7.627812308590531e-07, 2.7016636896519775e-07, 7.591345867972363e-08, 1.8882039831131855e-08, 6.461347548593791e-09, 5.169225787640778e-09, 4.7291286173112145e-09, 2.8045225869571647e-09, 9.93368827776587e-11, 3.738407450056794e-09, 2.242347297644185e-08, 5.63618320041525e-08, 9.454287791446423e-08, 1.1630420679175998e-07, 1.0932152533366832e-07, 7.390402146886183e-08, 2.9868511236445555e-08, 2.152625203499346e-09, 7.789024721578119e-09, 4.447370768285094e-08, 9.168721705976398e-08, 1.2327215881864552e-07, 1.2103213874988774e-07, 8.776368310067445e-08, 4.090986072926849e-08, 6.792643238846288e-09, 2.19668864711491e-09, 2.782163633755817e-08, 6.635125777562617e-08, 9.881463894052946e-08, 1.0351578526098045e-07, 7.680084200218828e-08, 4.002757625766943e-08, 9.199233876330439e-09, 3.19115890490225e-10, 1.5611848467412966e-08, 4.410871342345218e-08, 6.891131806253674e-08, 7.557326604980442e-08, 6.169749743287123e-08, 3.473783407151327e-08, 1.0090183910022984e-08, 7.828829157793151e-12, 8.407947966440945e-09, 2.8637561802230894e-08, 4.845193112634732e-08, 5.6149958709493355e-08, 4.833789301084396e-08, 2.9418708026228782e-08, 1.0113557978984861e-08, 3.087655437950417e-10, 4.241781579308114e-09, 1.8526256856559443e-08, 3.3786335472489205e-08, 4.519306841951586e-08, 3.881448924698381e-08, 2.4719796087206734e-08]

DEBUG chan1_val=[0.0, 5.3854044863015365e-09, 9.640159368285085e-09, 3.349331768132017e-08, 1.2421795158253949e-07, 3.6341984948339783e-07, 8.23077077942907e-07, 1.434576716733853e-06, 1.866701645782782e-06, 1.5438542213310247e-06, 3.749338479457899e-07, 7.453507980024111e-07, 1.1734790414117223e-05, 5.616548444848669e-05, 0.00018068424567010128, 0.0004680435489184571, 0.0010496614977870032, 0.0021151477032320854, 0.0039148671233587866, 0.006752328823753575, 0.010962257572673802, 0.016876927357780173, 0.024780600108025415, 0.034858301325087, 0.04715905003895009, 0.06154428582301777, 0.07769345569525589, 0.09512741328641873, 0.11321600322713167, 0.1312666841409905, 0.1485829112656172, 0.16453937332574944, 0.17864245366993378, 0.19057111985898165, 0.2001917864213171, 0.2075527594453615, 0.21285035078881917, 0.21639064100043984, 0.21853746985129235, 0.219666915917596, 0.22013033981703115, 0.22022056995088893, 0.220163368843537, 0.2201031226035194, 0.22011659330139172, 0.22021797264144286, 0.22038108739816417, 0.22055673058745498, 0.2206930376254265, 0.2207372669569089, 0.22070758290120332, 0.22058538776000827, 0.22041213092480963, 0.2202432845233554, 0.2201282920211581, 0.22009909252043325, 0.22015011804164322, 0.22021622167009153, 0.22016566190056386, 0.21978685315372312, 0.21879773896591448, 0.21685559767124882, 0.21358617058202842, 0.20861977210641258, 0.201637555300252, 0.19241945348148715, 0.18089024224208397, 0.1671489874054408, 0.1514856501802982, 0.13436490874145035, 0.11639353647841032, 0.09826009609865981, 0.08066416721162886, 0.06424851798518862, 0.049526982326139043, 0.03684895420515964, 0.02638152006067827, 0.018108002786643852, 0.011864160884747819, 0.007379407848727824, 0.00432669240423032, 0.0023686593438595836, 0.0011946184992459707, 0.000543853999802344, 0.00021612325412225832, 7.034810979818844e-05, 1.6143144968937013e-05, 1.4721664363097808e-06, 1.6130677541276803e-07, 1.3364463710935434e-06, 1.8251910984518694e-06, 1.4912391904610041e-06, 8.928723271066026e-07, 4.132461263639473e-07, 1.4927032080376913e-07, 4.30588437268206e-08, 1.2694930080039359e-08, 7.177480933041283e-09, 7.558659443929833e-09, 6.699986288340145e-09, 4.260144745690583e-09]

DEBUG chan2_val=[0.0, 3.641565425418014e-08, 2.1704963927836704e-08, 6.922057362411836e-09, 1.2047007023132046e-10, 5.340188635141553e-09, 1.9687948349177564e-08, 3.549742482142242e-08, 4.351839678999071e-08, 3.972930308789286e-08, 2.556572161499496e-08, 9.405954259029415e-09, 4.0264881028948915e-10, 4.245338640365165e-09, 1.934257373725204e-08, 3.769154249691701e-08, 4.910659829479462e-08, 4.704313410056672e-08, 3.25577183311115e-08, 1.3537276788669605e-08, 1.1948757556742458e-09, 3.1937990982054544e-09, 1.9612865752669968e-08, 4.1587587524942335e-08, 5.921265716191834e-08, 6.011841470022969e-08, 4.310321548306845e-08, 2.044273147120427e-08, 3.0300988390801734e-09, 1.872283197657291e-09, 1.846664132836599e-08, 4.4343113227461746e-08, 6.543763857457678e-08, 6.990710629098873e-08, 5.527016885013935e-08, 2.9392906460833003e-08, 6.9584352701171325e-09, 2.6919739111387274e-10, 1.2616675395223678e-08, 3.6218173074329675e-08, 5.7523544696087626e-08, 6.453111022487143e-08, 5.486669614767991e-08, 3.450003002160012e-08, 1.4576985916381805e-08, 2.87869110878829e-09, 5.959941221998039e-10, 3.364292991166032e-09, 6.148892223036055e-09, 7.611576000014078e-09, 7.121714024780303e-09, 1.0045074796507302e-08, 3.0689382708770264e-08, 1.0983173957770112e-07, 3.237270849963872e-07, 7.465631536417474e-07, 1.3368895672909508e-06, 1.7938644253757863e-06, 1.5623990088836515e-06, 4.4581471120703453e-07, 4.766312730569111e-07, 9.901338559772957e-06, 4.997441063509816e-05, 0.00016478702183477637, 0.00043343818144111866, 0.0009826315310130767, 0.0019966824121218844, 0.003720834638481129, 0.006454516741120437, 0.010531023490266127, 0.016284358578908496, 0.024004688346554705, 0.03388869922235809, 0.04599501965916217, 0.06020723914886733, 0.07622279515931596, 0.09356738968022632, 0.11162651435544899, 0.1297088761104789, 0.14711555742341745, 0.16321270433749818, 0.1774931049642609, 0.18961979884477909, 0.19944302446509216, 0.20699547166008497, 0.21246258865714163, 0.21614233181750966, 0.2183956163035948, 0.21959965610824125, 0.2201083442106375, 0.2202217297402528, 0.22017009293468717, 0.2201064060291775, 0.22011172291611983, 0.22020604142747316, 0.22036573245293614, 0.2205426745068929, 0.2206845121998875, 0.22074579106975312, 0.22071639968821308, 0.2205989989684551]

The real issue is the code that is generating x here [15]:

chan0_val = []
chan1_val = []
chan2_val = []
first_time = True
for p in period:
    t = np.arange(0, M*P*W)
    x = np.sin(t * p)
    X_psd = pfb_spectrometer(x, n_taps=M, n_chan=P, n_int=256, window_fn="hamming")
    if first_time:
        first_time = False
        print("DEBUG ")
    chan0_val.append(X_psd[0, 0])
    chan1_val.append(X_psd[0, 1])
    chan2_val.append(X_psd[0, 2])

The first pass through the loop, X_psd is all zeros. That's why the first element of chan_{0,1,2}_val = 0 causing db() to throw an exception. So, why was x = 0? x = np.sin(t * p) and the first value of t = 0.

My simple-minded fix:
x = np.sin(t * p) + 0.001 # the old "fudge factor" as it were

Now, it sails through just fine.