Rutherford Roof 30MHz Re-channelization Experiment¶

Key takaways:

  • Least squares optimization (with no added prior (other than gaussian sampling of ADC)) doesn't significantly improve on the Wiener filter
  • By salvaging 5% of the time field samples we can improve our worst quantization RMSE by a significant factor. The more we salvage the cleaner our signal gets.
  • Tentatively: the wider our FFTs the better they channelize. We can see this from the RMSE of the naive IPFBs are better for higher U. This is probably because the random bursts of noise at regular intervals cancel each other if there are many of them in the FFT as they are random, so the longer our FFTs the more our signal rises above the noise.
In [1]:
import numpy as np
from numpy.fft import fftfreq,fftshift
import matplotlib.pyplot as plt
import sys
sys.path.append("..") # for pfb tools imports etc

from pfb import forward_pfb,inverse_pfb,quantize_8_bit,quantize_8_bit_real,quantize_8_bit_spec_scaled_per_channel
import conjugate_gradient as cg
In [2]:
# shortcut variables
sqrt=np.sqrt

Global Variables¶

In [3]:
LFRAME=2048        # Length of a U1 frame
NCHAN=LFRAME//2+1  # Number of channels
NTAP=4             # Number of taps
NFRAMES=1600*32    # Number of frames of data to inspect
NFRAMES_OFFSET=100 # Number of frames to skip at the beginning
SR=10e6            # 10MHz Sample Rate
CF=35e6            # 35MHz Central Frequency 
FNAME="gqrx_20230601_175822_35000000_10000000_fc.raw" 
                   # File name/Relative path to data file ^
ADC_DELTA=0.04     # ADC quantization interval
COMPLEX_DELTA=0.25 # Quantization delta of channelized signal
IMSCALE_PLT=(-2.5,1)# Set scale of spectrogram plot colors
IMSCALE_RES=(-2.5,-1)# Set scale of residual spectrogram colors
CMAP_RES="terrain" # colormap of residual plts, default "gist_rainbow"
LEN_SIG=NFRAMES*LFRAME# length of entire signal

Waterfall plotting tool

In [4]:
from mpl_toolkits.axes_grid1 import make_axes_locatable
def waterfall(spec,sr,cf,title=None,imshow_scale=IMSCALE_PLT,
              cmap="gist_rainbow"):
    """Waterfall plot.
    
    Parameters
    ----------
    spec : 2darray
        The spectrum to plot
    sr : float
        Sample Rate (ADC)
    cf : float
        Central Frequency
    title : str
        The title of the plot. 
    imshow_scale : tuple
        Two tuple of floats. The first one sets the min 
        value to display, the second sets max value.
    cmap : str
        Name of the pyplot colormap.
    """
    if title is None:
        title=f"Central Freq {cf/1e6:.0f}MHz, Sample Rate {sr/1e6:.0f}MHz"
    nframes,nchan=spec.shape
    lframe=2*(nchan-1)
    fig,ax=plt.subplots(figsize=(6,6))
    log_spec = np.log(abs(spec)+1e-16) # add some noise for plotting purposes
    
    # Set the color scale by "modifying" spectrum (purely for visuals)
    log_spec[0,0]=imshow_scale[1] # Hack, set first pixel to cieling
    log_spec[0,1]=imshow_scale[0] # Hack, set second pixel to floor
    if (log_spec>imshow_scale[1]).any():
        print("WARNiNG: Spectrum has higher ceiling than is displayed!")
    log_spec[np.where(log_spec>imshow_scale[1])] = imshow_scale[1]
    log_spec[np.where(log_spec<imshow_scale[0])] = imshow_scale[0]
    
    freqs=fftshift(fftfreq(lframe,d=1/sr))[::2] + cf
    im = ax.imshow(log_spec,cmap=cmap)
    plt.title(title,fontsize=18)
    ax.set_xlabel("Frequency MHz")
    # plt.xticks(freqs)
    plt.axis('auto')
    xticks=np.arange(0,len(freqs),len(freqs)//6)
    ax.set_xticks(xticks)
    ax.set_xticklabels(["{:.1f}".format(i) for i in freqs[xticks]/1e6])

    ax.set_ylabel("Time (s)")
    yticks=np.arange(0,nframes,nframes//6)
    
    # colorbar stuff
    divider=make_axes_locatable(ax)
    cax = divider.append_axes('right',size='5%',pad=0.05)
    fig.colorbar(im,cax,orientation='vertical')
    
    ax.set_yticks(yticks)
    ax.set_yticklabels(["{:.2f}".format(i) for i in yticks * lframe/sr])
    
    plt.tight_layout()
    plt.show()

Load the signal¶

In [7]:
# Load from raw file
def load_sig(path:str, verbose=False):
    """Load the signal from binary file located at path"""
    sig=np.fromfile(path,dtype="float32") # load signal from file
    sig=(sig-sig.mean())/sig.std() # normalize
    sig=sig[LFRAME*NFRAMES_OFFSET:(NFRAMES+NFRAMES_OFFSET)*LFRAME] # trim signal
    sigadc=quantize_8_bit_real(sig,delta=ADC_DELTA) # ADC quantize signal
    if verbose is True:
        plt.figure()
        plt.plot(sig[:512],label="RAW")
        plt.plot(sigadc[:512],linewidth=1.5,label="ADC")
        plt.legend()
        plt.tight_layout()
        plt.show()
    return sigadc
SIG=load_sig(FNAME)
assert len(SIG)<=LEN_SIG, "not enough data for truncation: {len(sig)} < {LEN_SIG}"
SIG=SIG[:LEN_SIG] # trunkate signal

Channelize and Quantize¶

In [10]:
# Get the full channelized signal (for later)
spec_adc_full=forward_pfb(SIG,NCHAN,NTAP)
# Quantize the spectrum, 4+4bits=1byte (real+imag)
SPEC_Q,std_per_chan=quantize_8_bit_spec_scaled_per_channel(
        spec_adc_full, COMPLEX_DELTA)
del spec_adc_full # don't need this, takes up lots of memory
del std_per_chan # don't need this, cleaner code
In [11]:
def display_channelized_signal(sig,u=1):
    print("INFO: channelizing spectrum")
    spec_adc=forward_pfb(sig, (NCHAN-1)*u+1, NTAP)    
    print("INFO: Plotting channelized spectrum")
    waterfall(spec_adc/sqrt(u*LFRAME), # normalize spectrum
              SR, CF, title=f"Channelized ADC'd E-field U={u}")

    print("INFO: 4+4 bit quantizing spectrum")
    spec_adc_q,_=quantize_8_bit_spec_scaled_per_channel(
                    spec_adc, COMPLEX_DELTA)
    print("INFO: Plotting residuals: 'quantized spec' - 'truth'")
    waterfall((spec_adc_q-spec_adc)/spec_adc.std(axis=0),
                SR, CF, title=f"Channelized ADC'd E-field U={u}\
\n4+4 bit-quantized, Residuals\
\n(normalized by power in each channel)",
                imshow_scale=IMSCALE_RES, cmap=CMAP_RES)
    
    print("INFO: Garbage collection")
    del _
    del spec_adc_q
    del spec_adc
    return
In [12]:
display_channelized_signal(SIG,u=1)

INFO: channelizing spectrum
INFO: Plotting channelized spectrum
WARNiNG: Spectrum has higher ceiling than is displayed!

INFO: 4+4 bit quantizing spectrum
INFO: Plotting residuals: 'quantized spec' - 'truth'
WARNiNG: Spectrum has higher ceiling than is displayed!

INFO: Garbage collection
In [13]:
display_channelized_signal(SIG,u=2)

INFO: channelizing spectrum
INFO: Plotting channelized spectrum
WARNiNG: Spectrum has higher ceiling than is displayed!

INFO: 4+4 bit quantizing spectrum
INFO: Plotting residuals: 'quantized spec' - 'truth'
WARNiNG: Spectrum has higher ceiling than is displayed!

INFO: Garbage collection
In [14]:
display_channelized_signal(SIG,u=4)

INFO: channelizing spectrum
INFO: Plotting channelized spectrum
WARNiNG: Spectrum has higher ceiling than is displayed!

INFO: 4+4 bit quantizing spectrum
INFO: Plotting residuals: 'quantized spec' - 'truth'
WARNiNG: Spectrum has higher ceiling than is displayed!

INFO: Garbage collection
In [15]:
display_channelized_signal(SIG,u=8)

INFO: channelizing spectrum
INFO: Plotting channelized spectrum
WARNiNG: Spectrum has higher ceiling than is displayed!

INFO: 4+4 bit quantizing spectrum
INFO: Plotting residuals: 'quantized spec' - 'truth'
WARNiNG: Spectrum has higher ceiling than is displayed!

INFO: Garbage collection
In [16]:
display_channelized_signal(SIG,u=16)

INFO: channelizing spectrum
INFO: Plotting channelized spectrum
WARNiNG: Spectrum has higher ceiling than is displayed!

INFO: 4+4 bit quantizing spectrum
INFO: Plotting residuals: 'quantized spec' - 'truth'
WARNiNG: Spectrum has higher ceiling than is displayed!

INFO: Garbage collection
In [17]:
display_channelized_signal(SIG,u=32)

INFO: channelizing spectrum
INFO: Plotting channelized spectrum
WARNiNG: Spectrum has higher ceiling than is displayed!

INFO: 4+4 bit quantizing spectrum
INFO: Plotting residuals: 'quantized spec' - 'truth'
WARNiNG: Spectrum has higher ceiling than is displayed!

INFO: Garbage collection

Re-channelize naivly¶

Invert the PFB naively (without filter)¶
In [18]:
# def get_quantized_spec_and_ipfb(sig,u=2,wiener_thresh=0.1):
#     print("INFO: Channelizing signal")
#     spec_adc=forward_pfb(sig, (NCHAN-1)*u+1, NTAP)
#     print("INFO: 4+4 bit quantizing spectrum")
#     spec_adc_q,_=quantize_8_bit_spec_scaled_per_channel(
#                     spec_adc, COMPLEX_DELTA)
#     # garbage collction
#     del spec_adc
#     del _ 
#     print("INFO: INverting quantized spectrum")
#     sig_ipfb=inverse_pfb(spec_adc_q,wiener_thresh=wiener_thresh)
#     # drop imaginary componant
#     sig_ipfb,imag = sig_ipfb.real,sig_ipfb.imag 
#     assert sig_ipfb.std()/imag.std() > 1e10 # sanity check
#     print("DEBUG: Imaginary componant is negligible, cast to real")
#     del imag
#     return spec_adc_q,sig_ipfb
In [ ]:
In [19]:
# wrapper for inverse_pfb
def get_sig_ipfb(spec,wiener_thresh=0.0):
    print("INFO: Inverting spectrum")
    sig_ipfb=inverse_pfb(spec,wiener_thresh=wiener_thresh)
    # drop imaginary componant
    sig_ipfb,imag = sig_ipfb.real,sig_ipfb.imag 
    assert sig_ipfb.std()/imag.std() > 1e10 # sanity check
    print("DEBUG: Imaginary componant is negligible, cast to real")
    del imag
    return sig_ipfb # real signal
In [20]:
# Thin wrapper for forward pfb
def get_spec(sig, u):
    print("INFO: Channelizing signal")
    spec=forward_pfb(sig, (NCHAN-1)*u+1, NTAP)
    return spec

# get the quantized spectrum
def get_spec_q(sig, u):
    spec=get_spec(sig, u)
    print("INFO: 4+4 bit quantizing spectrum")
    spec_q,_=quantize_8_bit_spec_scaled_per_channel(
                    spec, COMPLEX_DELTA)
    del spec
    del _
    return spec_q
In [21]:
# Display spec then display residuals in waterfall plot
def display_spec_and_resid(spec, spec_true, u, title="Spec"):
    print("INFO: Display channelized spec")
    waterfall(spec/sqrt(u*LFRAME), SR, CF, title=title)
    print("INFO: Display residuals")
    waterfall((spec-spec_true)/spec_true.std(axis=0),
              SR, CF, title="Residuals, normalized by channel",
              imshow_scale=IMSCALE_RES, cmap=CMAP_RES)
    return
In [22]:
def rechannelize_ipfb(sig,sig_ipfb,u=2):
    spec_q = get_spec_q(sig, u)
    print(f"INFO: Re-channelize, upchan factor U={u}")
    spec_ipfb=get_spec(sig_ipfb, u)
    display_spec_and_resid(spec_ipfb, spec_q, u, title=f"Re-channelized IPFB'd signal U={u}")
    # Garbage collection
    del spec_ipfb
    del spec_q
    return
In [23]:
# Invert the signal without any filtering, FFT IPFB
sig_ipfb=get_sig_ipfb(SPEC_Q, wiener_thresh=0.0) # wiener_thresh=0 means no wiener filter
rechannelize_ipfb(SIG, sig_ipfb, u=1)

INFO: Inverting spectrum
DEBUG: Imaginary componant is negligible, cast to real
INFO: Channelizing signal
INFO: 4+4 bit quantizing spectrum
INFO: Re-channelize, upchan factor U=1
INFO: Channelizing signal
INFO: Display channelized spec
WARNiNG: Spectrum has higher ceiling than is displayed!

INFO: Display residuals
WARNiNG: Spectrum has higher ceiling than is displayed!
In [24]:
rechannelize_ipfb(SIG, sig_ipfb, u=2)

INFO: Channelizing signal
INFO: 4+4 bit quantizing spectrum
INFO: Re-channelize, upchan factor U=2
INFO: Channelizing signal
INFO: Display channelized spec
WARNiNG: Spectrum has higher ceiling than is displayed!

INFO: Display residuals
WARNiNG: Spectrum has higher ceiling than is displayed!
In [25]:
rechannelize_ipfb(SIG, sig_ipfb, u=4)

INFO: Channelizing signal
INFO: 4+4 bit quantizing spectrum
INFO: Re-channelize, upchan factor U=4
INFO: Channelizing signal
INFO: Display channelized spec
WARNiNG: Spectrum has higher ceiling than is displayed!

INFO: Display residuals
WARNiNG: Spectrum has higher ceiling than is displayed!
In [26]:
rechannelize_ipfb(SIG, sig_ipfb, u=8)

INFO: Channelizing signal
INFO: 4+4 bit quantizing spectrum
INFO: Re-channelize, upchan factor U=8
INFO: Channelizing signal
INFO: Display channelized spec
WARNiNG: Spectrum has higher ceiling than is displayed!

INFO: Display residuals
WARNiNG: Spectrum has higher ceiling than is displayed!
In [27]:
rechannelize_ipfb(SIG, sig_ipfb, u=16)

INFO: Channelizing signal
INFO: 4+4 bit quantizing spectrum
INFO: Re-channelize, upchan factor U=16
INFO: Channelizing signal
INFO: Display channelized spec
WARNiNG: Spectrum has higher ceiling than is displayed!

INFO: Display residuals
WARNiNG: Spectrum has higher ceiling than is displayed!
In [28]:
rechannelize_ipfb(SIG, sig_ipfb, u=32)

INFO: Channelizing signal
INFO: 4+4 bit quantizing spectrum
INFO: Re-channelize, upchan factor U=32
INFO: Channelizing signal
INFO: Display channelized spec
WARNiNG: Spectrum has higher ceiling than is displayed!

INFO: Display residuals
WARNiNG: Spectrum has higher ceiling than is displayed!
In [29]:
# Garbage collection
del sig_ipfb

re-channelize with Wiener filter¶

In [30]:
sig_wien=get_sig_ipfb(SPEC_Q, wiener_thresh=0.1)
rechannelize_ipfb(SIG, sig_wien, u=1)

INFO: Inverting spectrum
DEBUG: Imaginary componant is negligible, cast to real
INFO: Channelizing signal
INFO: 4+4 bit quantizing spectrum
INFO: Re-channelize, upchan factor U=1
INFO: Channelizing signal
INFO: Display channelized spec
WARNiNG: Spectrum has higher ceiling than is displayed!

INFO: Display residuals
WARNiNG: Spectrum has higher ceiling than is displayed!
In [31]:
rechannelize_ipfb(SIG, sig_wien, u=2)

INFO: Channelizing signal
INFO: 4+4 bit quantizing spectrum
INFO: Re-channelize, upchan factor U=2
INFO: Channelizing signal
INFO: Display channelized spec
WARNiNG: Spectrum has higher ceiling than is displayed!

INFO: Display residuals
WARNiNG: Spectrum has higher ceiling than is displayed!
In [32]:
rechannelize_ipfb(SIG, sig_wien, u=4)

INFO: Channelizing signal
INFO: 4+4 bit quantizing spectrum
INFO: Re-channelize, upchan factor U=4
INFO: Channelizing signal
INFO: Display channelized spec
WARNiNG: Spectrum has higher ceiling than is displayed!

INFO: Display residuals
WARNiNG: Spectrum has higher ceiling than is displayed!
In [33]:
rechannelize_ipfb(SIG, sig_wien, u=8)

INFO: Channelizing signal
INFO: 4+4 bit quantizing spectrum
INFO: Re-channelize, upchan factor U=8
INFO: Channelizing signal
INFO: Display channelized spec
WARNiNG: Spectrum has higher ceiling than is displayed!

INFO: Display residuals
WARNiNG: Spectrum has higher ceiling than is displayed!
In [34]:
rechannelize_ipfb(SIG, sig_wien, u=16)

INFO: Channelizing signal
INFO: 4+4 bit quantizing spectrum
INFO: Re-channelize, upchan factor U=16
INFO: Channelizing signal
INFO: Display channelized spec
WARNiNG: Spectrum has higher ceiling than is displayed!

INFO: Display residuals
WARNiNG: Spectrum has higher ceiling than is displayed!
In [35]:
rechannelize_ipfb(SIG, sig_wien, u=32)

INFO: Channelizing signal
INFO: 4+4 bit quantizing spectrum
INFO: Re-channelize, upchan factor U=32
INFO: Channelizing signal
INFO: Display channelized spec
WARNiNG: Spectrum has higher ceiling than is displayed!

INFO: Display residuals
WARNiNG: Spectrum has higher ceiling than is displayed!

Re-channelize with extra information¶

From IPFB Wiener, do CG-descent with no prior¶
In [36]:
# Construct B operator and u vector without a prior
B,u = cg.get_Bu_noprior(SPEC_Q,
                        delta=COMPLEX_DELTA,
                        ntap=NTAP)

# CG descent onto reconstructed timestream without prior
sig_cg_noprior=cg.conjugate_gradient_descent(B,u,x0=sig_wien,
                                             x_true=SIG,
                                             max_iter=5, 
                                             verbose=True,
                                             k=LEN_SIG//LFRAME,
                                             rmin=0)

# Garbage collection
del sig_cg_noprior 
del B
del u

INFO: Conjugate Gradient descent completed.
From IPFB Wiener, do CG-descent with 5% extra info¶
In [39]:
# Load 5% indices that make up your prior
# (with 5%, best result npersave=2, max_iter=1)
_,saved_idxs_5perc=cg.get_saved_idxs(n_per_save=2, prct=0.05,
                                   k=NFRAMES, lblock=LFRAME)

del _ # garbage collection

# Construct B operator and u vector using prior
B,u,chisq = cg.get_Bu_withprior(SIG, SPEC_Q,
                                delta=COMPLEX_DELTA, 
                                saved_idxs=saved_idxs_5perc,
                                ntap=NTAP)

# CG descend onto reconstructed timestream `ipfb_cg`
sig_cg_5perc=cg.conjugate_gradient_descent(B,u,x0=sig_wien,
                                           x_true=SIG,
                                           max_iter=1,
                                           verbose=True,
                                           k=NFRAMES,
                                           rmin=0)

# Garbage collection
del B
del u
del chisq # unused callable

INFO: Conjugate Gradient descent completed.
In [40]:
# Display channelized signal, U=1
for u in (1,2,4,8,16,32):
    spec_cg_5perc=get_spec(sig_cg_5perc, u)
    spec_adc = get_spec(SIG, u) # get true spectrum
    display_spec_and_resid(spec=spec_cg_5perc, spec_true=spec_adc, u=u, 
                           title=f"Spectrum CG-descent 5% extra data U={u}")
    del spec_cg_5perc # garbage collection
    del spec_adc

INFO: Channelizing signal
INFO: Channelizing signal
INFO: Display channelized spec
WARNiNG: Spectrum has higher ceiling than is displayed!

INFO: Display residuals
WARNiNG: Spectrum has higher ceiling than is displayed!

INFO: Channelizing signal
INFO: Channelizing signal
INFO: Display channelized spec
WARNiNG: Spectrum has higher ceiling than is displayed!

INFO: Display residuals
WARNiNG: Spectrum has higher ceiling than is displayed!

INFO: Channelizing signal
INFO: Channelizing signal
INFO: Display channelized spec
WARNiNG: Spectrum has higher ceiling than is displayed!

INFO: Display residuals
WARNiNG: Spectrum has higher ceiling than is displayed!

INFO: Channelizing signal
INFO: Channelizing signal
INFO: Display channelized spec
WARNiNG: Spectrum has higher ceiling than is displayed!

INFO: Display residuals
WARNiNG: Spectrum has higher ceiling than is displayed!

INFO: Channelizing signal
INFO: Channelizing signal
INFO: Display channelized spec
WARNiNG: Spectrum has higher ceiling than is displayed!

INFO: Display residuals
WARNiNG: Spectrum has higher ceiling than is displayed!

INFO: Channelizing signal
INFO: Channelizing signal
INFO: Display channelized spec
WARNiNG: Spectrum has higher ceiling than is displayed!

INFO: Display residuals
WARNiNG: Spectrum has higher ceiling than is displayed!
In [43]:
del sig_cg_5perc

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Input In [43], in <cell line: 1>()
----> 1 del sig_cg_5perc

NameError: name 'sig_cg_5perc' is not defined
From IPFB Wiener, do CG-descent with 10% extra info¶
In [45]:
# Load 5% indices that make up your prior
# (with 5%, best result npersave=1, max_iter=1 (confusingly means 1step))
_,saved_idxs_10perc=cg.get_saved_idxs(n_per_save=1, prct=0.10,
                                   k=NFRAMES, lblock=LFRAME)

# Construct B operator and u vector using prior
B,u,chisq = cg.get_Bu_withprior(SIG, SPEC_Q,
                                delta=COMPLEX_DELTA, 
                                saved_idxs=saved_idxs_10perc,
                                ntap=NTAP)

# CG descend onto reconstructed timestream `ipfb_cg`
sig_cg_10perc=cg.conjugate_gradient_descent(B,u,x0=sig_wien,
                                           x_true=SIG,
                                           max_iter=1,
                                           verbose=True,
                                           k=NFRAMES,
                                           rmin=0)

del B
del u
del chisq # unused callable

INFO: Conjugate Gradient descent completed.
In [46]:
# Display channelized signal
for u in (1,2,4,8,16,32):
    spec_cg_10perc=get_spec(sig_cg_10perc, u)
    spec_adc = get_spec(SIG, u) # get true spectrum
    display_spec_and_resid(spec=spec_cg_10perc, spec_true=spec_adc, u=u, 
                           title=f"Spectrum CG-descent 10% extra data U={u}")
    del spec_cg_10perc # garbage collection
    del spec_adc

INFO: Channelizing signal
INFO: Channelizing signal
INFO: Display channelized spec
WARNiNG: Spectrum has higher ceiling than is displayed!

INFO: Display residuals
WARNiNG: Spectrum has higher ceiling than is displayed!

INFO: Channelizing signal
INFO: Channelizing signal
INFO: Display channelized spec
WARNiNG: Spectrum has higher ceiling than is displayed!

INFO: Display residuals
WARNiNG: Spectrum has higher ceiling than is displayed!

INFO: Channelizing signal
INFO: Channelizing signal
INFO: Display channelized spec
WARNiNG: Spectrum has higher ceiling than is displayed!

INFO: Display residuals
WARNiNG: Spectrum has higher ceiling than is displayed!

INFO: Channelizing signal
INFO: Channelizing signal
INFO: Display channelized spec
WARNiNG: Spectrum has higher ceiling than is displayed!

INFO: Display residuals
WARNiNG: Spectrum has higher ceiling than is displayed!

INFO: Channelizing signal
INFO: Channelizing signal
INFO: Display channelized spec
WARNiNG: Spectrum has higher ceiling than is displayed!

INFO: Display residuals
WARNiNG: Spectrum has higher ceiling than is displayed!

INFO: Channelizing signal
INFO: Channelizing signal
INFO: Display channelized spec
WARNiNG: Spectrum has higher ceiling than is displayed!

INFO: Display residuals
WARNiNG: Spectrum has higher ceiling than is displayed!
In [51]:
# Garbage collection
del saved_idxs_10perc
del sig_cg_10perc

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Input In [51], in <cell line: 2>()
      1 # Garbage collection
----> 2 del saved_idxs_10perc
      3 del sig_cg_10perc

NameError: name 'saved_idxs_10perc' is not defined
From IPFB Wiener, do CG-descent with 20% extra info¶
In [48]:
# Load 5% indices that make up your prior
# (with 5%, best result npersave=2, max_iter=2 (confusingly means 1step))
_,saved_idxs_20perc=cg.get_saved_idxs(n_per_save=1, prct=0.20,
                                   k=NFRAMES, lblock=LFRAME)

# Construct B operator and u vector using prior
B,u,chisq = cg.get_Bu_withprior(SIG, SPEC_Q,
                                delta=COMPLEX_DELTA, 
                                saved_idxs=saved_idxs_20perc,
                                ntap=NTAP)

# CG descend onto reconstructed timestream `ipfb_cg`
sig_cg_20perc=cg.conjugate_gradient_descent(B,u,x0=sig_wien,
                                           x_true=SIG,
                                           max_iter=1,
                                           verbose=True,
                                           k=NFRAMES,
                                           rmin=0)

del B
del u
del chisq

INFO: Conjugate Gradient descent completed.
In [49]:
# Display channelized signal
for u in (1,2,4,8,16,32):
    spec_cg_20perc=get_spec(sig_cg_20perc, u)
    spec_adc = get_spec(SIG, u) # get true spectrum
    display_spec_and_resid(spec=spec_cg_20perc, spec_true=spec_adc, u=u, 
                           title=f"Spectrum CG-descent 20% extra data U={u}")
    del spec_cg_20perc # garbage collection
    del spec_adc

INFO: Channelizing signal
INFO: Channelizing signal
INFO: Display channelized spec
WARNiNG: Spectrum has higher ceiling than is displayed!

INFO: Display residuals
WARNiNG: Spectrum has higher ceiling than is displayed!

INFO: Channelizing signal
INFO: Channelizing signal
INFO: Display channelized spec
WARNiNG: Spectrum has higher ceiling than is displayed!

INFO: Display residuals
WARNiNG: Spectrum has higher ceiling than is displayed!

INFO: Channelizing signal
INFO: Channelizing signal
INFO: Display channelized spec
WARNiNG: Spectrum has higher ceiling than is displayed!

INFO: Display residuals
WARNiNG: Spectrum has higher ceiling than is displayed!

INFO: Channelizing signal
INFO: Channelizing signal
INFO: Display channelized spec
WARNiNG: Spectrum has higher ceiling than is displayed!

INFO: Display residuals
WARNiNG: Spectrum has higher ceiling than is displayed!

INFO: Channelizing signal
INFO: Channelizing signal
INFO: Display channelized spec
WARNiNG: Spectrum has higher ceiling than is displayed!

INFO: Display residuals
WARNiNG: Spectrum has higher ceiling than is displayed!

INFO: Channelizing signal
INFO: Channelizing signal
INFO: Display channelized spec
WARNiNG: Spectrum has higher ceiling than is displayed!

INFO: Display residuals
WARNiNG: Spectrum has higher ceiling than is displayed!
In [50]:
# Garbage collection
del saved_idxs_20perc
del sig_cg_20perc

Compare and Contrast Residuals¶

Below blocks no longer run because we've deleted the data from variables (otherwise RAM gets full and SWAP is used)

up-channelize by 2x, i.e. U=2¶

In [ ]:
# # IPFB
# waterfall((spec_ipfb_u2-spec_adc_u2_q)/spec_adc_u2_q.std(axis=0),
#           SR, CF, title=f"Residuals, IPFB'd naively\
# \nnormalized by channel",
#           imshow_scale=IMSCALE_RES,cmap=CMAP_RES)

# # Wiener filter
# waterfall((spec_wien_u2-spec_adc_u2_q)/spec_adc_u2_q.std(axis=0),
#           SR, CF, title=f"Residuals, Wiener Filtered\
# \nnormalized by channel",
#           imshow_scale=IMSCALE_RES,cmap=CMAP_RES)

# # 5% extra timestream data saved
# waterfall((spec_cg_5perc_u2-spec_adc_u2_q)/spec_adc_u2_q.std(axis=0),
#           SR, CF,
#           title=f"Residuals, CG 5% data saved\
# \nnormalized by channel",
#           imshow_scale=IMSCALE_RES,cmap=CMAP_RES)

# # 10% extra timestream data saved
# waterfall((spec_cg_10perc_u2-spec_adc_u2_q)/spec_adc_u2_q.std(axis=0),
#           SR, CF,
#           title=f"Residuals, CG 10% data saved\
# \nnormalized by channel",
#           imshow_scale=IMSCALE_RES,cmap=CMAP_RES)

# # 20% extra timestream data saved
# waterfall((spec_cg_20perc_u2-spec_adc_u2_q)/spec_adc_u2_q.std(axis=0),
#           SR, CF,
#           title=f"Residuals, CG 20% data saved\
# \nnormalized by channel",
#           imshow_scale=IMSCALE_RES,cmap=CMAP_RES)

up-channelize by 4x, i.e. U=4¶

In [ ]:
# # IPFB
# waterfall((spec_ipfb_u4-spec_adc_u4_q)/spec_adc_u4_q.std(axis=0),
#           SR, CF, title=f"Residuals, IPFB'd naively\
# \nnormalized by channel",
#           imshow_scale=IMSCALE_RES,cmap=CMAP_RES)

# # Wiener filter
# waterfall((spec_wien_u4-spec_adc_u4_q)/spec_adc_u4_q.std(axis=0),
#           SR, CF, title=f"Residuals, Wiener Filtered\
# \nnormalized by channel",
#           imshow_scale=IMSCALE_RES,cmap=CMAP_RES)

# # 5% extra timestream data saved
# waterfall((spec_cg_5perc_u4-spec_adc_u4_q)/spec_adc_u4_q.std(axis=0),
#           SR, CF,
#           title=f"Residuals, CG 5% data saved\
# \nnormalized by channel",
#           imshow_scale=IMSCALE_RES,cmap=CMAP_RES)

# # 10% extra timestream data saved
# waterfall((spec_cg_10perc_u4-spec_adc_u4_q)/spec_adc_u4_q.std(axis=0),
#           SR, CF,
#           title=f"Residuals, CG 10% data saved\
# \nnormalized by channel",
#           imshow_scale=IMSCALE_RES,cmap=CMAP_RES)

# # 20% extra timestream data saved
# waterfall((spec_cg_20perc_u4-spec_adc_u4_q)/spec_adc_u4_q.std(axis=0),
#           SR, CF,
#           title=f"Residuals, CG 20% data saved\
# \nnormalized by channel",
#           imshow_scale=IMSCALE_RES,cmap=CMAP_RES)
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]: