telegraphic / pfb_introduction

An interactive introduction to the polyphase filterbank technique for radio astronomy spectrometers
51 stars 16 forks source link

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

Closed texadactyl closed 4 years ago

texadactyl commented 4 years ago

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, Xpsd 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.

telegraphic commented 4 years ago

Agreed a fudge factor is nicer than a divide by zero error 👍