<jupyter_start><jupyter_text>Custom Signal InjectionThis tutorial shows how to add custom signal sources to synthetic voltage data. 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'
%load_ext autoreload
%autoreload 2

%matplotlib inline
import matplotlib.pyplot as plt

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

import sys
sys.path.insert(0, "../../")
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>We can create the backend elements up front, before the `Antenna` and signals.<jupyter_code>sample_rate = 3e9
num_taps = 8
num_branches = 1024

chan_bw = sample_rate / num_branches

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)<jupyter_output><empty_output><jupyter_text>Adding a custom signalLikewise, we can create the backend object before actually adding signals and noise to the `Antenna`, since data is only calculated when we call the `record` method. This way, we can actually use parameters specific to the background / recording in order to figure out how to construct the signal in the first place.<jupyter_code>antenna = stg.voltage.Antenna(sample_rate=sample_rate, 
                              fch1=6*u.GHz,
                              ascending=True,
                              num_pols=1)

rvb = stg.voltage.RawVoltageBackend(antenna,
                                    digitizer=digitizer,
                                    filterbank=filterbank,
                                    requantizer=requantizer,
                                    start_chan=0,
                                    num_chans=64,
                                    block_size=134217728,
                                    blocks_per_file=128,
                                    num_subblocks=32)<jupyter_output><empty_output><jupyter_text>We first add some Gaussian noise:<jupyter_code>antenna.x.add_noise(0, 1)

antenna.x.noise_std<jupyter_output><empty_output><jupyter_text>Custom signalsSimilarly, instead of the using `add_signal` method, we can add signal "sources", functions that accept `ts` as input and add the appropriate real voltages to the antenna data streams. Here's an example, utilizing the `scipy.signal` library. Alternatively, if you have `cusignal` installed, the injection will also be GPU accelerated.Note that if you have CuPy installed, arrays will be cast to CuPy arrays. Nevertheless, if custom signal source functions are written in NumPy, the data generation step won't be automatically GPU accelerated. So when possible, use CuPy to write functions. Keep in mind that CuPy is essentially a stand-in for NumPy, as is `cusignal` for `scipy.signal`. Of course, you can also easily manually convert between CuPy and NumPy via `cupy.asnumpy` and `cupy.asarray`.Generating a quadratic signal, scaled to amplitude 0.004, as well as a custom noise source:<jupyter_code>try:
    from cusignal import chirp
except ImportError:
    from scipy.signal import chirp

def my_signal(ts):
    return 0.004 * chirp(t=ts,
                         f0=6002.2e6,
                         t1=rvb.time_per_block,
                         f1=6002.3e6,
                         method='quadratic')

antenna.x.add_signal(my_signal)

def my_noise_source(ts):
    return np.random.normal(0, 2, len(ts))

antenna.x.add_signal(my_noise_source)<jupyter_output><empty_output><jupyter_text>However, after adding custom noise sources, the saved noise standard deviation may no longer be accurate.<jupyter_code>antenna.x.noise_std<jupyter_output><empty_output><jupyter_text>To fix this, we can run `DataStream.update_noise()` to estimate the noise statistics based on a few initial voltage samples. We can verify that the noise variance matches 1 + 2^2 = 5, as expected:<jupyter_code>antenna.x.update_noise()

antenna.x.noise_std**2<jupyter_output><empty_output><jupyter_text>Recording synthetic dataNow that we have all the pieces of the backend created, we can now create GUPPI RAW files using `RawVoltageBackend.record()`. There are two modes to record, `num_blocks` and `obs_length`. The former allows you to specify how many data blocks to record, as defined by the RAW format, and the latter computes the number of required blocks from the desired observation length (in seconds). The `header_dict` parameter allows for overwriting existing header values / adding additional keywords. For example, TurboSETI requires the `TELESCOP` parameter to be set to "GBT" to run on certain types of RAW files, so you can overwrite this parameter using `header_dict`.`raw_file_stem` is the path + prefix for your saved GUPPI RAW file.<jupyter_code>rvb.record(raw_file_stem='example_custom_signal',
           num_blocks=1, 
           length_mode='num_blocks',
           header_dict={'HELLO': 'test_value',
                        'TELESCOP': 'GBT'},
           verbose=False)<jupyter_output>Blocks:   0%|          | 0/1 [00:00<?, ?it/s]
  0%|          | 0/32 [00:00<?, ?it/s][A
Subblocks:   0%|          | 0/32 [00:00<?, ?it/s][A
Subblocks:   3%|▎         | 1/32 [00:02<01:27,  2.83s/it][A
Subblocks:   6%|▋         | 2/32 [00:04<01:13,  2.44s/it][A
Subblocks:   9%|▉         | 3/32 [00:05<01:01,  2.12s/it][A
Subblocks:  12%|█▎        | 4/32 [00:07<00:52,  1.89s/it][A
Subblocks:  16%|█▌        | 5/32 [00:08<00:46,  1.74s/it][A
Subblocks:  19%|█▉        | 6/32 [00:09<00:42,  1.63s/it][A
Subblocks:  22%|██▏       | 7/32 [00:11<00:38,  1.55s/it][A
Subblocks:  25%|██▌       | 8/32 [00:12<00:35,  1.49s/it][A
Subblocks:  28%|██▊       | 9/32 [00:13<00:33,  1.46s/it][A
Subblocks:  31%|███▏      | 10/32 [00:15<00:31,  1.43s/it][A
Subblocks:  34%|███▍      | 11/32 [00:16<00:29,  1.41s/it][A
Subblocks:  38%|███▊      | 12/32 [00:18<00:28,  1.40s/it][A
Subblocks:  41%|████      | 13/32 [00:19<00:26,  1.40s/it][A
Subblocks:  44%|████▍     | 14/32 [00:20<00:24,  1.39s/it][A
[...]<jupyter_text>Taking a lookIf `rawspec` is installed, you can unpack the data as such:<jupyter_code>!rawspec -f 1024 -t 16 -d . example_custom_signal<jupyter_output>working stem: example_custom_signal
opening file: example_custom_signal.0000.raw
opening file: example_custom_signal.0001.raw [No such file or directory]
output product 0: 64 spectra<jupyter_text>Using `blimpy` to open the resulting filterbank file and plotting the data:<jupyter_code>wf = bl.Waterfall('example_custom_signal.rawspec.0000.fil',
                  f_start=6002.0,
                  f_stop=6002.5)

plt.figure(figsize=(10, 6))
wf.plot_waterfall()
plt.show()<jupyter_output><empty_output><jupyter_text>Integrating over time:<jupyter_code>def db(x):
    """ Convert linear value to dB value """
    return 10*np.log10(x)

x = db(np.sum(wf.data[:, 0, :], axis=0))

plt.figure(figsize=(8, 6))
plt.plot(x)
plt.xlabel('Frequency bins')
plt.ylabel('Integrated Power (dB)')
plt.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.10494695129745622, 'fraction_numerical': 0.046657292598747284, 'mean_word_length': 3.2132471728594507, 'pattern_counts': {'":': 0, '<': 60, '<?xml version=': 0, '>': 43, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 13, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '11', 'dir': 'jupyter-structured-clean-dedup', 'id': '202006', 'license': 'MIT', 'n_tokens_mistral': 2844, 'n_tokens_neox': 2533, 'n_words': 745, 'path': 'jupyter-notebooks/voltage/05_custom_signals_estimate_noise.ipynb', 'repo_name': 'RocketRoss/setigen', 'size': '7823'}