Source code for NuRadioReco.detector.filterresponse

import numpy as np
from NuRadioReco.utilities import units
from scipy import interpolate as intp
import os
import glob


[docs]def get_filter_response_mini_circuits(frequencies, filter_name): """ Simple minicircuits filters. Check filter directory for available filters. """ # parse input file and convert to default units directory = os.path.dirname(os.path.abspath(__file__)) filename = os.path.join(directory, 'filter/{}_S2P'.format(filter_name), '{}_Minus40degC.S2P'.format(filter_name)) if not os.path.exists(filename): filename = os.path.join(directory, 'filter/{}_S2P'.format(filter_name), '{}__UNIT_1.s2p'.format(filter_name)) if not os.path.exists(filename): filename = glob.glob(os.path.join(directory, 'filter/{}_S2P'.format(filter_name), '{}*'.format(filter_name)))[0] data = np.loadtxt(filename, comments=['#', '!']) ff = data.T[0] * units.MHz S21 = 10**(data.T[3] / 20.) * np.exp(1j * data.T[4] * units.deg) get_S21 = intp.interp1d(ff, S21) response = np.zeros_like(frequencies, dtype=complex) mask = (frequencies > ff.min()) & (frequencies < ff.max()) response[mask] = get_S21(frequencies[mask]) return response
[docs]def get_filter_response_mini_circuits2(frequencies, filter_name): """ Simple minicircuits filters. Check filter directory for available filters. Parse input file and convert to default units """ directory = os.path.dirname(os.path.abspath(__file__)) ff, insertion_loss, return_loss, ff2, group_delay = np.loadtxt( os.path.join(directory, 'filter/{}.dat'.format(filter_name)), unpack=True) ff *= units.MHz ff2 *= units.MHz group_delay *= units.ns insertion_loss = 10**(-insertion_loss / 20.) get_insertion_loss = intp.interp1d(ff, insertion_loss) get_group_delay = intp.interp1d(ff2, group_delay) df = 1 * units.MHz fff2 = np.linspace(ff2.min(), ff2.max(), int(np.ceil((ff2.max() - ff2.min()) / (1 * units.MHz)))) phase2 = -2 * np.pi * np.cumsum(get_group_delay(fff2) * df) get_phase = intp.interp1d(fff2, phase2) response = np.zeros_like(frequencies, dtype=complex) mask = (frequencies > max(ff.min(), ff2.min())) & (frequencies < min(ff.max(), ff2.max())) response[mask] = get_insertion_loss(frequencies[mask]) * np.exp(1j * get_phase(frequencies[mask])) return response
[docs]def get_filter_response(frequencies, filter_name): """ Get measured filter responses. """ directory = os.path.dirname(os.path.abspath(__file__)) if filter_name == 'NTU+cheb': ff, mag, phase = np.loadtxt(os.path.join(directory, 'filter/NTU+cheb_filter_mag_phase.txt'), unpack=True) get_phase = intp.interp1d(ff, np.unwrap(phase)) get_insertion_loss = intp.interp1d(ff, mag) ff2 = ff else: # parse input file and convert to default units ff2, phase, t = np.loadtxt(os.path.join( directory, 'filter/measurement/{}_PHASE.CSV'.format(filter_name)), unpack=True, skiprows=3, delimiter=',') ff, insertion_loss, t = np.loadtxt(os.path.join( directory, 'filter/measurement/{}_LINMAG.CSV'.format(filter_name)), unpack=True, skiprows=3, delimiter=',') ff *= units.Hz ff2 *= units.Hz phase *= units.deg get_phase = intp.interp1d(ff2, np.unwrap(phase)) get_insertion_loss = intp.interp1d(ff, insertion_loss) response = np.zeros_like(frequencies, dtype=complex) mask = (frequencies > max(ff.min(), ff2.min())) & (frequencies < min(ff.max(), ff2.max())) response[mask] = get_insertion_loss(frequencies[mask]) * np.exp(1j * get_phase(frequencies[mask])) return response
if __name__ == "__main__": from matplotlib import pyplot as plt freqs = np.linspace(10 * units.MHz, 500 * units.MHz, 100) delta_f = freqs[1] - freqs[0] fig, (ax, ax2) = plt.subplots(1, 2) for name in ["SHP-48+", "SHP-20+", "SHP-25+", "SXHP-48+", "ZFHP-1R2-S+", "ZX75HP-44+"]: response1 = get_filter_response_mini_circuits(freqs, name) ax.plot(freqs / units.MHz, np.abs(response1), label=name) grp_delay = -1. / 2 / np.pi * np.diff(np.unwrap(np.angle(response1))) / delta_f # ax2.plot(ff / units.MHz, np.unwrap(np.angle(response1)) / units.deg) ax2.plot((freqs[1:] + freqs[:-1]) * 0.5 / units.MHz, grp_delay / units.ns, label=name) ax.legend() ax2.legend() fig.tight_layout() plt.show() freqs = np.linspace(200 * units.MHz, 800 * units.MHz, 100) fig, (ax, ax2) = plt.subplots(1, 2) for name in ["ZX75LP-470+", "ZFLP-450+", "SBLP-467+"]: response1 = get_filter_response_mini_circuits(freqs, name) ax.plot(freqs / units.MHz, np.abs(response1), label=name) grp_delay = -1. / 2 / np.pi * np.diff(np.unwrap(np.angle(response1))) / delta_f # ax2.plot(ff / units.MHz, np.unwrap(np.angle(response1)) / units.deg) ax2.plot((freqs[1:] + freqs[:-1]) * 0.5 / units.MHz, grp_delay / units.ns, label=name) ax.legend() ax2.legend() fig.tight_layout() plt.show() response1 = get_filter_response_mini_circuits(freqs, "SHP-100+") response2 = get_filter_response(freqs, "SHP-100+") fig, (ax, ax2) = plt.subplots(1, 2) phase_1 = np.unwrap(np.angle(response1)) phase_1 = phase_1 / np.min(phase_1) phase_2 = np.unwrap(np.angle(response2)) phase_2 = phase_2 / np.min(phase_2) ax.plot(freqs / units.MHz, phase_1, label='spec sheet') ax.plot(freqs / units.MHz, phase_2, label='measurement') ax2.plot(freqs[1:] / units.MHz, np.diff(np.unwrap(np.angle(response1))), label='spec sheet') ax2.plot(freqs[1:] / units.MHz, np.diff(np.unwrap(np.angle(response2))), label='measurement') ax.legend() plt.show()