0
# Thresholding and Denoising
1

2
Signal thresholding functions for noise reduction, feature extraction, and signal enhancement with various thresholding methods for wavelet coefficient processing.
3

4
## Capabilities
5

6
### Basic Thresholding
7

8
Fundamental thresholding operations for coefficient processing.
9

10
```python { .api }
11
def threshold(data, value, mode: str = 'soft', substitute: float = 0):
12
    """
13
    Apply thresholding to data.
14
    
15
    Parameters:
16
    - data: Input array to threshold
17
    - value: Threshold value (scalar or array-like for per-element thresholds)
18
    - mode: Thresholding mode ('soft', 'hard', 'garrote', 'greater', 'less')  
19
    - substitute: Substitute value for thresholded elements (default: 0)
20
    
21
    Returns:
22
    Thresholded array of same shape as input
23
    """
24
```
25

26
#### Usage Examples
27

28
```python
29
import pywt
30
import numpy as np
31
import matplotlib.pyplot as plt
32

33
# Create test signal with noise
34
np.random.seed(42)
35
t = np.linspace(0, 1, 1000)
36
clean_signal = np.sin(2 * np.pi * 5 * t) + 0.5 * np.sin(2 * np.pi * 20 * t)
37
noise = 0.5 * np.random.randn(len(t))
38
noisy_signal = clean_signal + noise
39

40
print(f"Signal length: {len(noisy_signal)}")
41
print(f"Noise std: {np.std(noise):.3f}")
42

43
# Wavelet decomposition
44
coeffs = pywt.wavedec(noisy_signal, 'db8', level=6)
45
print(f"Decomposition levels: {len(coeffs) - 1}")
46

47
# Different thresholding modes
48
threshold_value = 0.1
49
modes = ['soft', 'hard', 'garrote', 'greater', 'less']
50

51
# Apply different thresholding modes to detail coefficients
52
thresholded_results = {}
53

54
for mode in modes:
55
    coeffs_thresh = coeffs.copy()
56
    
57
    # Apply thresholding to all detail coefficients (skip approximation)
58
    for i in range(1, len(coeffs_thresh)):
59
        coeffs_thresh[i] = pywt.threshold(coeffs_thresh[i], threshold_value, mode=mode)
60
    
61
    # Reconstruct signal
62
    reconstructed = pywt.waverec(coeffs_thresh, 'db8')
63
    thresholded_results[mode] = reconstructed
64

65
# Visualization
66
fig, axes = plt.subplots(3, 2, figsize=(15, 12))
67
axes = axes.flatten()
68

69
# Original signals
70
axes[0].plot(t, clean_signal, 'b-', label='Clean signal', linewidth=2)
71
axes[0].plot(t, noisy_signal, 'r-', alpha=0.7, label='Noisy signal')
72
axes[0].set_title('Original Signals')
73
axes[0].legend()
74
axes[0].grid(True)
75

76
# Thresholded results
77
for i, (mode, result) in enumerate(thresholded_results.items()):
78
    axes[i+1].plot(t, clean_signal, 'b-', alpha=0.7, label='Clean')
79
    axes[i+1].plot(t, result, 'g-', linewidth=2, label=f'{mode.capitalize()} threshold')
80
    axes[i+1].set_title(f'{mode.capitalize()} Thresholding')
81
    axes[i+1].legend()
82
    axes[i+1].grid(True)
83

84
plt.tight_layout()
85
plt.show()
86

87
# Compare thresholding effectiveness
88
def calculate_mse(signal1, signal2):
89
    """Calculate Mean Squared Error."""
90
    return np.mean((signal1 - signal2)**2)
91

92
def calculate_snr(clean, noisy):
93
    """Calculate Signal-to-Noise Ratio in dB."""
94
    signal_power = np.mean(clean**2)
95
    noise_power = np.mean((noisy - clean)**2)
96
    return 10 * np.log10(signal_power / noise_power)
97

98
print("\nThresholding Performance:")
99
print(f"Original SNR: {calculate_snr(clean_signal, noisy_signal):.2f} dB")
100
print()
101

102
for mode, result in thresholded_results.items():
103
    mse = calculate_mse(clean_signal, result)
104
    snr = calculate_snr(clean_signal, result)
105
    print(f"{mode.capitalize():>8}: MSE = {mse:.6f}, SNR = {snr:.2f} dB")
106

107
# Demonstrate coefficient-specific thresholding
108
coeffs_detailed = pywt.wavedec(noisy_signal, 'db8', level=6)
109

110
# Analyze coefficient statistics
111
print("\nCoefficient Statistics:")
112
for i, coeff in enumerate(coeffs_detailed):
113
    coeff_type = "Approximation" if i == 0 else f"Detail {len(coeffs_detailed)-i}"
114
    print(f"{coeff_type:>15}: std = {np.std(coeff):.4f}, max = {np.max(np.abs(coeff)):.4f}")
115

116
# Adaptive thresholding based on coefficient statistics
117
coeffs_adaptive = coeffs_detailed.copy()
118
for i in range(1, len(coeffs_adaptive)):  # Skip approximation
119
    # Use different threshold for each level
120
    sigma = np.std(coeffs_adaptive[i])
121
    adaptive_threshold = 2 * sigma  # 2-sigma rule
122
    coeffs_adaptive[i] = pywt.threshold(coeffs_adaptive[i], adaptive_threshold, mode='soft')
123

124
adaptive_result = pywt.waverec(coeffs_adaptive, 'db8')
125
adaptive_snr = calculate_snr(clean_signal, adaptive_result)
126
print(f"\nAdaptive thresholding SNR: {adaptive_snr:.2f} dB")
127
```
128

129
### Advanced Thresholding
130

131
More sophisticated thresholding methods.
132

133
```python { .api }
134
def threshold_firm(data, value_low, value_high):
135
    """
136
    Firm thresholding (intermediate between soft and hard).
137
    
138
    Parameters:
139
    - data: Input array to threshold
140
    - value_low: Lower threshold value
141
    - value_high: Upper threshold value (must be >= value_low)
142
    
143
    Returns:
144
    Firm thresholded array
145
    """
146
```
147

148
#### Usage Examples
149

150
```python
151
import pywt
152
import numpy as np
153
import matplotlib.pyplot as plt
154

155
# Create test signal with mixed noise characteristics
156
t = np.linspace(0, 2, 2000)
157
signal = (np.sin(2 * np.pi * 3 * t) +           # Low frequency component
158
          0.7 * np.sin(2 * np.pi * 15 * t) +    # Medium frequency component
159
          0.4 * np.sin(2 * np.pi * 50 * t))     # High frequency component
160

161
# Add different types of noise
162
gaussian_noise = 0.3 * np.random.randn(len(t))
163
impulse_noise = np.zeros_like(t)
164
impulse_indices = np.random.choice(len(t), size=50, replace=False)
165
impulse_noise[impulse_indices] = 2 * np.random.randn(50)
166

167
mixed_noisy = signal + gaussian_noise + impulse_noise
168

169
# Wavelet decomposition
170
coeffs = pywt.wavedec(mixed_noisy, 'db6', level=7)
171

172
# Compare soft, hard, and firm thresholding
173
threshold_low = 0.05
174
threshold_high = 0.15
175

176
# Soft thresholding
177
coeffs_soft = coeffs.copy()
178
for i in range(1, len(coeffs_soft)):
179
    coeffs_soft[i] = pywt.threshold(coeffs_soft[i], threshold_high, mode='soft')
180
reconstructed_soft = pywt.waverec(coeffs_soft, 'db6')
181

182
# Hard thresholding  
183
coeffs_hard = coeffs.copy()
184
for i in range(1, len(coeffs_hard)):
185
    coeffs_hard[i] = pywt.threshold(coeffs_hard[i], threshold_high, mode='hard')
186
reconstructed_hard = pywt.waverec(coeffs_hard, 'db6')
187

188
# Firm thresholding
189
coeffs_firm = coeffs.copy()
190
for i in range(1, len(coeffs_firm)):
191
    coeffs_firm[i] = pywt.threshold_firm(coeffs_firm[i], threshold_low, threshold_high)
192
reconstructed_firm = pywt.waverec(coeffs_firm, 'db6')
193

194
# Visualization
195
fig, axes = plt.subplots(2, 2, figsize=(15, 10))
196

197
axes[0, 0].plot(t, signal, 'b-', label='Clean signal', linewidth=2)
198
axes[0, 0].plot(t, mixed_noisy, 'r-', alpha=0.7, label='Noisy signal')
199
axes[0, 0].set_title('Original vs Noisy Signal')
200
axes[0, 0].legend()
201
axes[0, 0].grid(True)
202

203
axes[0, 1].plot(t, signal, 'b-', alpha=0.7, label='Clean')
204
axes[0, 1].plot(t, reconstructed_soft, 'g-', linewidth=2, label='Soft threshold')
205
axes[0, 1].set_title('Soft Thresholding')
206
axes[0, 1].legend()
207
axes[0, 1].grid(True)
208

209
axes[1, 0].plot(t, signal, 'b-', alpha=0.7, label='Clean')
210
axes[1, 0].plot(t, reconstructed_hard, 'r-', linewidth=2, label='Hard threshold')
211
axes[1, 0].set_title('Hard Thresholding')
212
axes[1, 0].legend()
213
axes[1, 0].grid(True)
214

215
axes[1, 1].plot(t, signal, 'b-', alpha=0.7, label='Clean')
216
axes[1, 1].plot(t, reconstructed_firm, 'm-', linewidth=2, label='Firm threshold')
217
axes[1, 1].set_title('Firm Thresholding')
218
axes[1, 1].legend()
219
axes[1, 1].grid(True)
220

221
plt.tight_layout()
222
plt.show()
223

224
# Performance comparison
225
def calculate_metrics(clean, denoised):
226
    """Calculate denoising metrics."""
227
    mse = np.mean((clean - denoised)**2)
228
    snr = 10 * np.log10(np.mean(clean**2) / np.mean((clean - denoised)**2))
229
    
230
    # Peak Signal-to-Noise Ratio
231
    max_val = np.max(clean)
232
    psnr = 20 * np.log10(max_val / np.sqrt(mse))
233
    
234
    return mse, snr, psnr
235

236
print("Thresholding Method Comparison:")
237
print(f"{'Method':<15} {'MSE':<10} {'SNR (dB)':<10} {'PSNR (dB)':<10}")
238
print("-" * 50)
239

240
original_mse, original_snr, original_psnr = calculate_metrics(signal, mixed_noisy)
241
print(f"{'Original':<15} {original_mse:<10.6f} {original_snr:<10.2f} {original_psnr:<10.2f}")
242

243
soft_mse, soft_snr, soft_psnr = calculate_metrics(signal, reconstructed_soft)
244
print(f"{'Soft':<15} {soft_mse:<10.6f} {soft_snr:<10.2f} {soft_psnr:<10.2f}")
245

246
hard_mse, hard_snr, hard_psnr = calculate_metrics(signal, reconstructed_hard)
247
print(f"{'Hard':<15} {hard_mse:<10.6f} {hard_snr:<10.2f} {hard_psnr:<10.2f}")
248

249
firm_mse, firm_snr, firm_psnr = calculate_metrics(signal, reconstructed_firm)
250
print(f"{'Firm':<15} {firm_mse:<10.6f} {firm_snr:<10.2f} {firm_psnr:<10.2f}")
251

252
# Demonstrate threshold selection strategies
253
def sure_threshold_estimate(coeffs, sigma=None):
254
    """Estimate threshold using Stein's Unbiased Risk Estimate (SURE)."""
255
    if sigma is None:
256
        # Estimate noise standard deviation from finest detail coefficients
257
        sigma = np.median(np.abs(coeffs[-1])) / 0.6745
258
    
259
    n = len(coeffs[-1])
260
    threshold = sigma * np.sqrt(2 * np.log(n))
261
    return threshold
262

263
def bayes_threshold_estimate(coeffs):
264
    """Simple Bayesian threshold estimate."""
265
    # Estimate based on coefficient statistics
266
    all_details = np.concatenate([coeff for coeff in coeffs[1:]])
267
    return np.std(all_details)
268

269
# Apply different threshold selection methods
270
sure_thresh = sure_threshold_estimate(coeffs)
271
bayes_thresh = bayes_threshold_estimate(coeffs)
272

273
print(f"\nThreshold Selection:")
274
print(f"SURE threshold: {sure_thresh:.4f}")
275
print(f"Bayesian threshold: {bayes_thresh:.4f}")
276
print(f"Manual threshold: {threshold_high:.4f}")
277

278
# Test SURE threshold
279
coeffs_sure = coeffs.copy()
280
for i in range(1, len(coeffs_sure)):
281
    coeffs_sure[i] = pywt.threshold(coeffs_sure[i], sure_thresh, mode='soft')
282
reconstructed_sure = pywt.waverec(coeffs_sure, 'db6')
283

284
sure_mse, sure_snr, sure_psnr = calculate_metrics(signal, reconstructed_sure)
285
print(f"{'SURE':<15} {sure_mse:<10.6f} {sure_snr:<10.2f} {sure_psnr:<10.2f}")
286
```
287

288
### 2D Image Thresholding
289

290
Thresholding applications for image denoising and processing.
291

292
#### Usage Examples
293

294
```python
295
import pywt
296
import numpy as np
297
import matplotlib.pyplot as plt
298

299
# Create test image
300
size = 256
301
x, y = np.mgrid[0:size, 0:size]
302
image = np.zeros((size, size))
303

304
# Add geometric shapes
305
image[64:192, 64:192] = 1.0  # Square
306
center_x, center_y = 128, 128
307
radius = 40
308
circle_mask = (x - center_x)**2 + (y - center_y)**2 <= radius**2
309
image[circle_mask] = 0.5
310

311
# Add texture and noise
312
texture = 0.1 * np.sin(2 * np.pi * x / 16) * np.cos(2 * np.pi * y / 16)
313
noise = 0.3 * np.random.randn(size, size)  
314
noisy_image = image + texture + noise
315

316
print(f"Image shape: {noisy_image.shape}")
317
print(f"Pixel value range: [{np.min(noisy_image):.3f}, {np.max(noisy_image):.3f}]")
318

319
# 2D wavelet decomposition
320
coeffs_2d = pywt.wavedec2(noisy_image, 'db4', level=4)
321
print(f"2D decomposition levels: {len(coeffs_2d) - 1}")
322

323
# Apply thresholding to 2D coefficients
324
def threshold_2d_coeffs(coeffs, threshold_value, mode='soft'):
325
    """Apply thresholding to 2D wavelet coefficients."""
326
    coeffs_thresh = [coeffs[0]]  # Keep approximation unchanged
327
    
328
    for level_coeffs in coeffs[1:]:
329
        cH, cV, cD = level_coeffs
330
        cH_thresh = pywt.threshold(cH, threshold_value, mode=mode)
331
        cV_thresh = pywt.threshold(cV, threshold_value, mode=mode)
332
        cD_thresh = pywt.threshold(cD, threshold_value, mode=mode)
333
        coeffs_thresh.append((cH_thresh, cV_thresh, cD_thresh))
334
    
335
    return coeffs_thresh
336

337
# Different threshold values
338
thresholds = [0.05, 0.1, 0.2, 0.3]
339
denoised_images = {}
340

341
for thresh in thresholds:
342
    coeffs_thresh = threshold_2d_coeffs(coeffs_2d, thresh, mode='soft')
343
    denoised = pywt.waverec2(coeffs_thresh, 'db4')
344
    denoised_images[thresh] = denoised
345

346
# Visualization
347
fig, axes = plt.subplots(2, 3, figsize=(18, 12))
348

349
axes[0, 0].imshow(image, cmap='gray')
350
axes[0, 0].set_title('Original Clean Image')
351
axes[0, 0].axis('off')
352

353
axes[0, 1].imshow(noisy_image, cmap='gray')
354
axes[0, 1].set_title('Noisy Image')
355
axes[0, 1].axis('off')
356

357
# Show wavelet decomposition
358
cA = coeffs_2d[0]
359
axes[0, 2].imshow(cA, cmap='gray')
360
axes[0, 2].set_title('Approximation (Level 4)')
361
axes[0, 2].axis('off')
362

363
# Show denoised results
364
for i, (thresh, denoised) in enumerate(denoised_images.items()):
365
    if i < 3:
366
        axes[1, i].imshow(denoised, cmap='gray')
367
        axes[1, i].set_title(f'Denoised (threshold = {thresh})')
368
        axes[1, i].axis('off')
369

370
plt.tight_layout()
371
plt.show()
372

373
# Quantitative evaluation
374
def image_quality_metrics(original, denoised):
375
    """Calculate image quality metrics."""
376
    mse = np.mean((original - denoised)**2)
377
    
378
    # Peak Signal-to-Noise Ratio
379
    max_val = np.max(original)
380
    psnr = 20 * np.log10(max_val / np.sqrt(mse)) if mse > 0 else float('inf')
381
    
382
    # Structural Similarity Index (simplified)
383
    mu1, mu2 = np.mean(original), np.mean(denoised)
384
    sigma1, sigma2 = np.std(original), np.std(denoised)
385
    sigma12 = np.mean((original - mu1) * (denoised - mu2))
386
    
387
    c1, c2 = 0.01**2, 0.03**2  # Constants for stability
388
    ssim = ((2 * mu1 * mu2 + c1) * (2 * sigma12 + c2)) / ((mu1**2 + mu2**2 + c1) * (sigma1**2 + sigma2**2 + c2))
389
    
390
    return mse, psnr, ssim
391

392
print("\n2D Image Denoising Results:")
393
print(f"{'Threshold':<10} {'MSE':<10} {'PSNR (dB)':<12} {'SSIM':<10}")
394
print("-" * 50)
395

396
for thresh, denoised in denoised_images.items():
397
    mse, psnr, ssim = image_quality_metrics(image, denoised)
398
    print(f"{thresh:<10} {mse:<10.6f} {psnr:<12.2f} {ssim:<10.4f}")
399

400
# Adaptive 2D thresholding
401
def adaptive_2d_threshold(coeffs):
402
    """Apply adaptive thresholding based on coefficient statistics at each level."""
403
    coeffs_adaptive = [coeffs[0]]  # Keep approximation
404
    
405
    for level, (cH, cV, cD) in enumerate(coeffs[1:]):
406
        # Calculate adaptive thresholds for each orientation
407
        thresh_H = 3 * np.std(cH)
408
        thresh_V = 3 * np.std(cV) 
409
        thresh_D = 3 * np.std(cD)
410
        
411
        cH_thresh = pywt.threshold(cH, thresh_H, mode='soft')
412
        cV_thresh = pywt.threshold(cV, thresh_V, mode='soft')
413
        cD_thresh = pywt.threshold(cD, thresh_D, mode='soft')
414
        
415
        coeffs_adaptive.append((cH_thresh, cV_thresh, cD_thresh))
416
        
417
        print(f"Level {len(coeffs)-1-level}: thresholds H={thresh_H:.4f}, V={thresh_V:.4f}, D={thresh_D:.4f}")
418
    
419
    return coeffs_adaptive
420

421
coeffs_adaptive = adaptive_2d_threshold(coeffs_2d)
422
adaptive_denoised = pywt.waverec2(coeffs_adaptive, 'db4')
423

424
adaptive_mse, adaptive_psnr, adaptive_ssim = image_quality_metrics(image, adaptive_denoised)
425
print(f"{'Adaptive':<10} {adaptive_mse:<10.6f} {adaptive_psnr:<12.2f} {adaptive_ssim:<10.4f}")
426

427
# Show adaptive result
428
plt.figure(figsize=(15, 5))
429
plt.subplot(1, 3, 1)
430
plt.imshow(image, cmap='gray')
431
plt.title('Original')
432
plt.axis('off')
433

434
plt.subplot(1, 3, 2)
435
plt.imshow(noisy_image, cmap='gray')
436
plt.title('Noisy')
437
plt.axis('off')
438

439
plt.subplot(1, 3, 3)
440
plt.imshow(adaptive_denoised, cmap='gray')
441
plt.title('Adaptive Denoised')
442
plt.axis('off')
443

444
plt.tight_layout()
445
plt.show()
446
```
447

448
## Types
449

450
```python { .api }
451
# Thresholding modes
452
ThresholdMode = Literal['soft', 'hard', 'garrote', 'greater', 'less']
453

454
# Threshold value types
455
ThresholdValue = Union[float, np.ndarray]  # Scalar or array for per-element thresholds
456

457
# Common threshold selection methods
458
ThresholdMethod = Literal['sure', 'bayes', 'minimax', 'manual']
459
```