Question about using DPLMS to get optimum filter for HPGe detector data

[zql2021]

Hello, Thank you very much for developing and maintaining this project. I often refer to it during digital signal processing.

I am analyzing data from HPGe detector (preamlifier output), after applying a trapezoidal filter, the energy resolution was not satisfactory. Therefore, I try the DPLMS method to synthesize an optimal filter, mainly referencing
src/pygama/pargen/dplms_ge_dict.py and https://doi.org/10.1140/epjc/s10052-023-11299-z

However, I have encountered some issues. Below are my processing steps and results. I would greatly appreciate any suggestions for checks or debugging.

1. Waveform basic info

Sampling rate: 100 MHz
Wave length: 60 µs per waveform
Baseline length: ~10 µs
Rise time: 100–300 ns
Here is an example of the experimental waveform

2. Data selection & baseline subtraction

Selected Cs-137 full-energy-peak waveforms. For each waveform, baseline is calculated and subtracted:

baseline_mean = np.mean(elec_waves[:, :900], axis=1)     
wave_subbase = elec_waves - baseline_mean[:, None]  

3. Noise matrix calculation

I took the first 9 µs baseline region to compute the noise matrix:

length = 200
bls_data = wave_subbase[:, :900]   
def noise_matrix(bls: np.array, length: int) -> np.array:
    nev, size = bls.shape
    nmat = np.matmul(bls.T, bls, dtype=float) / nev
    kernel = np.identity(size - length + 1)
    nmat = convolve2d(nmat, kernel, boundary="symm", mode="valid") / (size - length + 1)
    return nmat
nmat = noise_matrix(bls_data, length)

4. Signal alignment & reference matrix

Aligned waveforms using trigger position, then extracted ±3 µs around the trigger for wfs (size = 600). Here are the aligned waveforms.

Computed reference waveform and reference matrix（Flat-top and pile-up constraints are temporarily not considered.）

ref = np.mean(wfs, axis=0)   
ref /= np.max(ref)
flo = size // 2 
fhi = size // 2 + length 
ref_segment = ref[flo:fhi]
rmat = np.outer(ref_segment, ref_segment)  # (length, length)

and the ref waveform:

5. Filter synthesis

def filter_synthesis(ref, nmat, rmat, za, pmat, fmat, length, size, flip=True):
    flo = size // 2 
    fhi = size // 2 + length
    ref_window = ref[flo:fhi]
    mat = nmat + rmat + za * np.ones([length, length]) + pmat + fmat
    x = np.linalg.solve(mat, ref_window).astype(np.float32)
    y = convolve(ref, np.flip(x), mode="valid")
    maxy = np.max(y)
    x /= maxy
    y /= maxy
    refy = ref[(size // 2) - (len(y) // 2):(size // 2) + (len(y) // 2)]
    if flip:
        return np.flip(x), y, refy
    else:
        return x, y, refy

only the noise matrix (nmat) and the reference matrix (rmat) are considered for testing.

coeffs = {
    "nm": 0.5,   
    "za": 0., 
    "pl": 0,  
    "ft": 0.0    
}

length = nmat.shape[0]
size   = len(ref)
x, y, refy = filter_synthesis(
    ref,
    coeffs["nm"] * nmat,
    rmat,
    coeffs["za"],
    pmat * coeffs["pl"],
    fmat * coeffs["ft"],
    length,
    size,
    flip=True
)

6. Visualization

plt.figure(figsize=(10,4))
plt.subplot(1,2,1)
plt.plot(x, label="Filter coefficients")
plt.xlabel("Tap index")
plt.ylabel("Coefficient value")
plt.legend()
plt.title("DPLMS Filter Coefficients")

plt.subplot(1,2,2)
plt.plot(y, label="Filter output")
plt.plot(refy, label="Reference")
plt.xlabel("Sample index")
plt.ylabel("Amplitude")
plt.legend()
plt.title("Filtered output vs Reference")
plt.tight_layout()
plt.show()

Final results show that y (filtered output) and refy (reference) differ significantly in the rising edge:

Question:

What could be the possible reasons for the difference between y and refy in the rising edge of the waveform?
Could you suggest some debugging ideas or parameter adjustments that might help? Thank you in advance for your help!

[valerioda]

Thanks a lot for starting this discussion! From what I can see, the procedure looks mostly correct.

I'm a bit uncertain about the indices (flo and fhi) used to calculate ref_segment from ref, and later in the filter calculation.
In our implementation, they are defined as:

flo = size // 2 - length // 2
fhi = size // 2 + length // 2

I’d suggest using the functions already defined in the routine whenever possible, to ensure consistency.

A general comment about the filter: a flat top will be needed to properly calculate the energy.

For any further discussion, feel free to reach out at valerio.dandrea@lngs.infn.it