<jupyter_start><jupyter_text>Efficient RAW file generation by subsamplingThis tutorial discusses generating high resolution synthetic data with smaller volumes by using lower sample rates. If you have access to a GPU, it is highly recommended to install CuPy, which performs the equivalent NumPy array operations on the GPU (https://docs.cupy.dev/en/stable/install.html). This is not necessary to run raw voltage generation, but will highly accelerate the pipeline. Once you have CuPy installed, to enable GPU acceleration you must set `SETIGEN_ENABLE_GPU` to '1' in the shell or in Python via `os.environ`. It can also be useful to set `CUDA_VISIBLE_DEVICES` to specify which GPUs to use.<jupyter_code># !pip install cupy-cuda110
import os
os.environ['SETIGEN_ENABLE_GPU'] = '1'
os.environ['CUDA_VISIBLE_DEVICES'] = '0'
%matplotlib inline
import matplotlib.pyplot as plt

import numpy as np
from astropy import units as u
import blimpy as bl

import setigen as stg
# Sometimes it can be necessary to re-run this command for plots to show automatically
%matplotlib inline<jupyter_output><empty_output><jupyter_text>Especially when using a high sampling rate like 3 GHz, `setigen`'s raw voltage pipeline can produce a *lot* of data, most of which is empty. Using the standard synthesis method, you can throw out information after it's generated, but that still requires generating an unnecessarily large amount of data beforehand. How can we cut down on unnecessary computation, but maintain time and frequency resolutions, not to mention signal processing conversions?It turns out that as long as we're careful, we can reduce the sampling rate and the number of PFB branches and maintain the same overall structure of the data:<jupyter_code>subsample_factor = 128

sample_rate = 3e9 // subsample_factor
num_taps = 8
num_branches = 1024 // subsample_factor

print(f'Max "num_chans" is {num_branches // 2}.')<jupyter_output>Max "num_chans" is 4.<jupyter_text>Note that the number of coarse channels to record `num_chans` has a reduced maximal value, since we are subsampling. An additional effect when reducing the data volume by a signficiant factor is that if you keep the `block_size` the same, you may end up writing more data when recording in `num_blocks` mode. This might be the desired behavior, but it's important to keep in mind. Since `record()` will generate enough data to fill out full blocks, if you additionally want to reduce the amount of data recorded, be sure to accordingly reduce `block_size`. So to test this out, let's add the same kinds of signals as we did in `05_raw_file_gen_snr.ipynb`, so that we can directly compare synthetic data products.Setting up backend elements, then adding noise and a few signals:<jupyter_code>fftlength = 1024
int_factor = 4
num_blocks = 1

digitizer = stg.voltage.RealQuantizer(target_fwhm=32,
                                      num_bits=8)

filterbank = stg.voltage.PolyphaseFilterbank(num_taps=num_taps, 
                                             num_branches=num_branches)

requantizer = stg.voltage.ComplexQuantizer(target_fwhm=32,
                                           num_bits=8)

antenna = stg.voltage.Antenna(sample_rate=sample_rate, 
                              fch1=6*u.GHz,
                              ascending=True,
                              num_pols=2)

block_size = stg.voltage.get_block_size(num_antennas=1,
                                        tchans_per_block=128,
                                        num_bits=8,
                                        num_pols=2,
                                        num_branches=num_branches,
                                        num_chans=num_branches//2,
                                        fftlength=fftlength,
                                        int_factor=int_factor)
rvb = stg.voltage.RawVoltageBackend(antenna,
                                    digitizer=digitizer,
                                    filterbank=filterbank,
                                    requantizer=requantizer,
                                    start_chan=0,
                                    num_chans=num_branches//2,
                                    block_size=block_size,
                                    blocks_per_file=128,
                                    num_subblocks=32)

# Add noise
for stream in antenna.streams:
    stream.add_noise(v_mean=0, 
                     v_std=1)
    
# Add signals
signal_level = stg.voltage.get_level(10, 
                                     rvb,
                                     num_blocks=num_blocks,
                                     length_mode='num_blocks',
                                     fftlength=fftlength)

for f_start in np.linspace(6003.1e6, 6003.9e6, 9):
    leakage_factor = stg.voltage.get_leakage_factor(f_start, rvb, fftlength)
    for stream in antenna.streams:
        level = stream.get_total_noise_std() * leakage_factor * signal_level
        stream.add_constant_signal(f_start=f_start, 
                                   drift_rate=0*u.Hz/u.s, 
                                   level=level)
    
rvb.record(output_file_stem='example_subsample',
           num_blocks=num_blocks, 
           length_mode='num_blocks',
           header_dict={'HELLO': 'test_value',
                        'TELESCOP': 'GBT'},
           verbose=False)

!rawspec -f $fftlength -t $int_factor -d . example_subsample<jupyter_output>Blocks:   0%|          | 0/1 [00:00<?, ?it/s]
  0%|          | 0/64 [00:00<?, ?it/s][A
Subblocks:   0%|          | 0/64 [00:00<?, ?it/s][A
Subblocks:   2%|▏         | 1/64 [00:01<01:31,  1.45s/it][A
Subblocks:  12%|█▎        | 8/64 [00:01<00:57,  1.02s/it][A
Subblocks:  39%|███▉      | 25/64 [00:01<00:27,  1.39it/s][A
Subblocks:  66%|██████▌   | 42/64 [00:01<00:11,  1.98it/s][A
Subblocks:  92%|█████████▏| 59/64 [00:01<00:01,  2.82it/s][A
Blocks: 100%|██████████| 1/1 [00:02<00:00,  2.02s/it]     [A<jupyter_text>Plotting the results:<jupyter_code>wf = bl.Waterfall('example_subsample.rawspec.0000.fil',
                  f_start=6003.0,
                  f_stop=6004.0)

plt.figure(figsize=(10, 6))
wf.plot_waterfall()
plt.show()

frame = stg.Frame(wf)

spectrum = stg.integrate(frame, normalize=True)

plt.figure(figsize=(8, 6))
plt.plot(spectrum)
plt.xlabel('Frequency bins')
plt.ylabel('SNR')
plt.show()<jupyter_output><empty_output><jupyter_text>This looks very similar to the plot in `05_raw_file_gen_snr.ipynb`, as expected! We can also check the frame resolution:<jupyter_code>frame.get_params()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.0881057268722467, 'fraction_numerical': 0.031900349384778974, 'mean_word_length': 2.338742393509128, 'pattern_counts': {'":': 0, '<': 28, '<?xml version=': 0, '>': 19, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 26, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '5', 'dir': 'jupyter-structured-clean-dedup', 'id': '63083', 'license': 'MIT', 'n_tokens_mistral': 2154, 'n_tokens_neox': 1945, 'n_words': 562, 'path': 'jupyter-notebooks/voltage/07_efficient_gen_by_subsampling.ipynb', 'repo_name': 'bbrzycki/setigen', 'size': '6583'}