Search code examples
pythonimage-thresholdingpywavelets

Stein's Unbiased Risk Estimate (SURE) threshold calculation (Wavelets)


I want to perform wavelaet threshold processing on 1d signal using pywavelets. I wrote this code to find argmin of risk like it mentioned here (https://cdn.intechopen.com/pdfs/34936/InTech-Wavelet_denoising.pdf) to calculate SURE threshold, because I didn't find complete solution in python packages:

min_risk = np.inf
sorted_coeffs = np.sort(np.abs(np.concatenate(coeffs)))
for i in range(len(sorted_coeffs)):
    res = (n-i)*sigma**2 + i*(sigma**2+sorted_coeffs[i]**2)
    risk = np.sum([sorted_coeffs[j]**2 for j in range(i, n)]) - res
    if risk < min_risk:
        min_risk = risk
        T_s = sorted_coeffs[i]

But actually it performs too long, if I have lots of coefficients (about 0.5 million). Is it possible to calculate this threshold more effectively?


Solution

  • As far as I understood, there was a solution in pyyawt package: https://pyyawt.readthedocs.io/_modules/pyyawt/denoising.html#thselect