0
# Multiresolution Analysis
1
2
Multiresolution analysis (MRA) providing direct access to approximation components at each decomposition level for signal analysis and reconstruction with multiple transform backends.
3
4
## Capabilities
5
6
### 1D Multiresolution Analysis
7
8
Decomposition into approximation components at different scales for 1D signals.
9
10
```python { .api }
11
def mra(data, wavelet, level: int = None, axis: int = -1,
12
transform: str = 'swt', mode: str = 'periodization'):
13
"""
14
1D multiresolution analysis.
15
16
Parameters:
17
- data: Input 1D signal
18
- wavelet: Wavelet specification
19
- level: Number of decomposition levels (default: maximum possible)
20
- axis: Axis along which to perform MRA (default: -1)
21
- transform: Transform type ('swt' for stationary WT, 'dwt' for discrete WT)
22
- mode: Signal extension mode (default: 'periodization')
23
24
Returns:
25
List of approximation components [A1, A2, ..., An] where:
26
- A1: Finest scale approximation
27
- An: Coarsest scale approximation
28
"""
29
30
def imra(mra_coeffs):
31
"""
32
Inverse 1D MRA via summation.
33
34
Parameters:
35
- mra_coeffs: List of approximation components from mra()
36
37
Returns:
38
Reconstructed signal (sum of all approximation levels)
39
"""
40
```
41
42
#### Usage Examples
43
44
```python
45
import pywt
46
import numpy as np
47
import matplotlib.pyplot as plt
48
49
# Create test signal with multiple scales
50
np.random.seed(42)
51
t = np.linspace(0, 4, 2048)
52
signal = (2 * np.sin(2 * np.pi * 1 * t) + # Slow trend
53
1.5 * np.sin(2 * np.pi * 8 * t) + # Medium frequency
54
1 * np.sin(2 * np.pi * 32 * t) + # Fast oscillation
55
0.5 * np.sin(2 * np.pi * 100 * t) + # High frequency
56
0.3 * np.random.randn(len(t))) # Noise
57
58
print(f"Signal length: {len(signal)}")
59
60
# Perform multiresolution analysis
61
level = 8
62
mra_coeffs = pywt.mra(signal, 'db8', level=level, transform='swt')
63
print(f"MRA decomposition levels: {len(mra_coeffs)}")
64
print(f"Each approximation shape: {mra_coeffs[0].shape}")
65
66
# Verify perfect reconstruction
67
reconstructed = pywt.imra(mra_coeffs)
68
reconstruction_error = np.max(np.abs(signal - reconstructed))
69
print(f"MRA reconstruction error: {reconstruction_error:.2e}")
70
71
# Analyze frequency content at each scale
72
def analyze_frequency_content(approx_component, sampling_rate=1.0, title=""):
73
"""Analyze dominant frequency in approximation component."""
74
# Simple frequency analysis using FFT
75
fft_vals = np.fft.fft(approx_component)
76
freqs = np.fft.fftfreq(len(approx_component), 1/sampling_rate)
77
78
# Find dominant frequency (excluding DC component)
79
magnitude = np.abs(fft_vals[1:len(fft_vals)//2])
80
freq_bins = freqs[1:len(freqs)//2]
81
dominant_idx = np.argmax(magnitude)
82
dominant_freq = freq_bins[dominant_idx]
83
84
return dominant_freq, magnitude, freq_bins
85
86
sampling_rate = len(t) / (t[-1] - t[0]) # Samples per second
87
print(f"\nFrequency analysis (sampling rate: {sampling_rate:.1f} Hz):")
88
89
dominant_frequencies = []
90
for i, approx in enumerate(mra_coeffs):
91
dom_freq, _, _ = analyze_frequency_content(approx, sampling_rate)
92
dominant_frequencies.append(dom_freq)
93
scale_factor = 2**(i+1)
94
expected_max_freq = sampling_rate / (2 * scale_factor)
95
print(f"Level {i+1}: dominant freq = {dom_freq:.2f} Hz, max expected = {expected_max_freq:.2f} Hz")
96
97
# Visualization of multiresolution components
98
fig, axes = plt.subplots(min(6, len(mra_coeffs) + 1), 1, figsize=(15, 12))
99
100
# Original signal
101
axes[0].plot(t, signal, 'k-', linewidth=1)
102
axes[0].set_title('Original Signal')
103
axes[0].set_ylabel('Amplitude')
104
axes[0].grid(True, alpha=0.3)
105
106
# Plot first 5 approximation levels
107
for i in range(min(5, len(mra_coeffs))):
108
axes[i+1].plot(t, mra_coeffs[i], 'b-', linewidth=1.5)
109
axes[i+1].set_title(f'Approximation Level {i+1} (Scale 2^{i+1})')
110
axes[i+1].set_ylabel('Amplitude')
111
axes[i+1].grid(True, alpha=0.3)
112
113
axes[-1].set_xlabel('Time')
114
plt.tight_layout()
115
plt.show()
116
117
# Compare different transform backends
118
transforms = ['swt', 'dwt']
119
mra_comparison = {}
120
121
for transform_type in transforms:
122
try:
123
if transform_type == 'dwt':
124
# DWT-based MRA may have different requirements
125
mra_dwt = pywt.mra(signal, 'db8', level=6, transform='dwt')
126
mra_comparison[transform_type] = mra_dwt
127
else:
128
mra_comparison[transform_type] = mra_coeffs
129
130
print(f"\n{transform_type.upper()}-based MRA:")
131
mra_result = mra_comparison[transform_type]
132
print(f" Levels: {len(mra_result)}")
133
print(f" Component shapes: {[comp.shape for comp in mra_result[:3]]}...")
134
135
# Reconstruction test
136
reconstructed_comp = pywt.imra(mra_result)
137
if len(reconstructed_comp) == len(signal):
138
error = np.max(np.abs(signal - reconstructed_comp))
139
print(f" Reconstruction error: {error:.2e}")
140
else:
141
print(f" Reconstruction shape: {reconstructed_comp.shape} (vs {signal.shape})")
142
143
except Exception as e:
144
print(f"Error with {transform_type}: {e}")
145
146
# Adaptive signal analysis using MRA
147
def adaptive_denoising_mra(signal, wavelet, noise_level_estimate=None):
148
"""Adaptive denoising using MRA with scale-dependent thresholding."""
149
150
# Perform MRA
151
mra_components = pywt.mra(signal, wavelet, transform='swt')
152
153
if noise_level_estimate is None:
154
# Estimate noise from finest scale approximation
155
finest_approx = mra_components[0]
156
noise_level_estimate = np.std(finest_approx - np.median(finest_approx))
157
158
print(f"Estimated noise level: {noise_level_estimate:.4f}")
159
160
# Apply scale-dependent denoising
161
denoised_components = []
162
for i, component in enumerate(mra_components):
163
scale = 2**(i+1)
164
# Threshold decreases with scale (coarser scales need less denoising)
165
threshold = noise_level_estimate / np.sqrt(scale)
166
167
# Apply soft thresholding
168
denoised = pywt.threshold(component, threshold, mode='soft')
169
denoised_components.append(denoised)
170
171
print(f"Level {i+1}: threshold = {threshold:.4f}")
172
173
# Reconstruct denoised signal
174
return pywt.imra(denoised_components)
175
176
# Test adaptive denoising
177
noisy_signal = signal + 0.5 * np.random.randn(len(signal))
178
denoised_mra = adaptive_denoising_mra(noisy_signal, 'db8', noise_level_estimate=0.5)
179
180
# Compare with standard wavelet denoising
181
coeffs_std = pywt.wavedec(noisy_signal, 'db8', level=8)
182
coeffs_thresh = [coeffs_std[0]] # Keep approximation
183
for detail in coeffs_std[1:]:
184
threshold = 0.1 * np.std(detail)
185
coeffs_thresh.append(pywt.threshold(detail, threshold, mode='soft'))
186
denoised_std = pywt.waverec(coeffs_thresh, 'db8')
187
188
# Performance comparison
189
def calculate_snr(clean, noisy):
190
return 10 * np.log10(np.var(clean) / np.var(clean - noisy))
191
192
original_snr = calculate_snr(signal, noisy_signal)
193
mra_snr = calculate_snr(signal, denoised_mra)
194
std_snr = calculate_snr(signal, denoised_std)
195
196
print(f"\nDenoising Performance:")
197
print(f"Original SNR: {original_snr:.2f} dB")
198
print(f"MRA denoising SNR: {mra_snr:.2f} dB")
199
print(f"Standard denoising SNR: {std_snr:.2f} dB")
200
201
# Plot denoising comparison
202
plt.figure(figsize=(15, 10))
203
204
plt.subplot(2, 2, 1)
205
plt.plot(t, signal, 'b-', label='Clean signal', linewidth=2)
206
plt.plot(t, noisy_signal, 'r-', alpha=0.7, label='Noisy signal')
207
plt.title('Original vs Noisy Signal')
208
plt.legend()
209
plt.grid(True)
210
211
plt.subplot(2, 2, 2)
212
plt.plot(t, signal, 'b-', alpha=0.7, label='Clean')
213
plt.plot(t, denoised_mra, 'g-', linewidth=2, label='MRA denoised')
214
plt.title(f'MRA Denoising (SNR: {mra_snr:.2f} dB)')
215
plt.legend()
216
plt.grid(True)
217
218
plt.subplot(2, 2, 3)
219
plt.plot(t, signal, 'b-', alpha=0.7, label='Clean')
220
plt.plot(t, denoised_std, 'm-', linewidth=2, label='Standard denoised')
221
plt.title(f'Standard Denoising (SNR: {std_snr:.2f} dB)')
222
plt.legend()
223
plt.grid(True)
224
225
plt.subplot(2, 2, 4)
226
plt.plot(t, signal - denoised_mra, 'g-', label='MRA error')
227
plt.plot(t, signal - denoised_std, 'm-', alpha=0.7, label='Standard error')
228
plt.title('Denoising Errors')
229
plt.legend()
230
plt.grid(True)
231
232
plt.tight_layout()
233
plt.show()
234
```
235
236
### 2D Multiresolution Analysis
237
238
Multiresolution decomposition for 2D images providing scale-space analysis.
239
240
```python { .api }
241
def mra2(data, wavelet, level: int = None, axes=(-2, -1),
242
transform: str = 'swt2', mode: str = 'periodization'):
243
"""
244
2D multiresolution analysis.
245
246
Parameters:
247
- data: Input 2D array (image)
248
- wavelet: Wavelet specification
249
- level: Number of decomposition levels (default: maximum possible)
250
- axes: Pair of axes for 2D transform (default: last two axes)
251
- transform: Transform type ('swt2' for 2D stationary WT)
252
- mode: Signal extension mode (default: 'periodization')
253
254
Returns:
255
List of 2D approximation components at different scales
256
"""
257
258
def imra2(mra_coeffs):
259
"""
260
Inverse 2D MRA via summation.
261
262
Parameters:
263
- mra_coeffs: List of 2D approximation components from mra2()
264
265
Returns:
266
Reconstructed 2D array (sum of all approximation levels)
267
"""
268
```
269
270
#### Usage Examples
271
272
```python
273
import pywt
274
import numpy as np
275
import matplotlib.pyplot as plt
276
277
# Create test image with multiple scales
278
size = 256
279
x, y = np.mgrid[0:size, 0:size]
280
281
# Multi-scale image with different frequency content
282
image = (0.5 + # DC component
283
0.4 * np.sin(2 * np.pi * x / 64) * np.cos(2 * np.pi * y / 64) + # Low frequency
284
0.3 * np.sin(2 * np.pi * x / 16) * np.cos(2 * np.pi * y / 16) + # Medium frequency
285
0.2 * np.sin(2 * np.pi * x / 4) * np.cos(2 * np.pi * y / 4)) # High frequency
286
287
# Add texture and noise
288
texture = 0.1 * np.random.randn(size, size)
289
image = image + texture
290
291
print(f"Image shape: {image.shape}")
292
print(f"Image value range: [{np.min(image):.3f}, {np.max(image):.3f}]")
293
294
# 2D Multiresolution analysis
295
level = 5
296
mra2_coeffs = pywt.mra2(image, 'db4', level=level, transform='swt2')
297
print(f"2D MRA levels: {len(mra2_coeffs)}")
298
print(f"Approximation shapes: {[approx.shape for approx in mra2_coeffs]}")
299
300
# Perfect reconstruction
301
reconstructed_2d = pywt.imra2(mra2_coeffs)
302
reconstruction_error_2d = np.max(np.abs(image - reconstructed_2d))
303
print(f"2D MRA reconstruction error: {reconstruction_error_2d:.2e}")
304
305
# Visualize 2D multiresolution components
306
fig, axes = plt.subplots(2, 3, figsize=(18, 12))
307
axes = axes.flatten()
308
309
# Original image
310
axes[0].imshow(image, cmap='gray')
311
axes[0].set_title('Original Image')
312
axes[0].axis('off')
313
314
# Show first 5 approximation levels
315
for i in range(min(5, len(mra2_coeffs))):
316
axes[i+1].imshow(mra2_coeffs[i], cmap='gray')
317
scale = 2**(i+1)
318
axes[i+1].set_title(f'Approximation Level {i+1} (Scale 2^{i+1})')
319
axes[i+1].axis('off')
320
321
plt.tight_layout()
322
plt.show()
323
324
# Scale-space analysis
325
def analyze_image_scales(mra_components):
326
"""Analyze image content at different scales."""
327
analysis = {}
328
329
for i, component in enumerate(mra_components):
330
scale = 2**(i+1)
331
332
# Calculate statistics
333
mean_val = np.mean(component)
334
std_val = np.std(component)
335
energy = np.sum(component**2)
336
337
# Edge content (simple gradient measure)
338
grad_x = np.diff(component, axis=1)
339
grad_y = np.diff(component, axis=0)
340
edge_strength = np.mean(np.sqrt(grad_x[:-1,:]**2 + grad_y[:,:-1]**2))
341
342
analysis[i+1] = {
343
'scale': scale,
344
'mean': mean_val,
345
'std': std_val,
346
'energy': energy,
347
'edge_strength': edge_strength,
348
'shape': component.shape
349
}
350
351
return analysis
352
353
scale_analysis = analyze_image_scales(mra2_coeffs)
354
355
print("\n2D Scale-Space Analysis:")
356
print(f"{'Level':<6} {'Scale':<6} {'Mean':<8} {'Std':<8} {'Energy':<10} {'Edges':<8}")
357
print("-" * 55)
358
359
for level, stats in scale_analysis.items():
360
print(f"{level:<6} {stats['scale']:<6} {stats['mean']:<8.4f} {stats['std']:<8.4f} "
361
f"{stats['energy']:<10.2f} {stats['edge_strength']:<8.4f}")
362
363
# Multi-scale image enhancement
364
def multiscale_enhancement(image, wavelet, enhancement_factors=None):
365
"""Enhance image using multi-scale processing."""
366
367
mra_components = pywt.mra2(image, wavelet, transform='swt2')
368
369
if enhancement_factors is None:
370
# Default: enhance fine details more than coarse
371
enhancement_factors = [1.5, 1.3, 1.1, 1.0, 0.9]
372
373
# Ensure we have enough factors
374
enhancement_factors = enhancement_factors[:len(mra_components)]
375
enhancement_factors += [1.0] * (len(mra_components) - len(enhancement_factors))
376
377
# Apply enhancement
378
enhanced_components = []
379
for i, (component, factor) in enumerate(zip(mra_components, enhancement_factors)):
380
if factor != 1.0:
381
# Enhance by scaling around the mean
382
mean_val = np.mean(component)
383
enhanced = mean_val + factor * (component - mean_val)
384
else:
385
enhanced = component
386
387
enhanced_components.append(enhanced)
388
print(f"Level {i+1}: enhancement factor = {factor}")
389
390
return pywt.imra2(enhanced_components)
391
392
# Apply enhancement
393
enhanced_image = multiscale_enhancement(image, 'db4', [2.0, 1.5, 1.2, 1.0, 0.8])
394
395
# Compare original and enhanced
396
plt.figure(figsize=(15, 5))
397
398
plt.subplot(1, 3, 1)
399
plt.imshow(image, cmap='gray')
400
plt.title('Original Image')
401
plt.axis('off')
402
403
plt.subplot(1, 3, 2)
404
plt.imshow(enhanced_image, cmap='gray')
405
plt.title('Enhanced Image')
406
plt.axis('off')
407
408
plt.subplot(1, 3, 3)
409
plt.imshow(enhanced_image - image, cmap='RdBu_r')
410
plt.title('Enhancement Difference')
411
plt.colorbar()
412
plt.axis('off')
413
414
plt.tight_layout()
415
plt.show()
416
417
# Calculate enhancement metrics
418
original_edge_strength = np.mean(np.abs(np.gradient(image)[0]) + np.abs(np.gradient(image)[1]))
419
enhanced_edge_strength = np.mean(np.abs(np.gradient(enhanced_image)[0]) + np.abs(np.gradient(enhanced_image)[1]))
420
421
print(f"\nEnhancement Results:")
422
print(f"Original edge strength: {original_edge_strength:.4f}")
423
print(f"Enhanced edge strength: {enhanced_edge_strength:.4f}")
424
print(f"Enhancement ratio: {enhanced_edge_strength / original_edge_strength:.2f}")
425
```
426
427
### nD Multiresolution Analysis
428
429
Multiresolution analysis for n-dimensional data.
430
431
```python { .api }
432
def mran(data, wavelet, level: int = None, axes=None,
433
transform: str = 'swtn', mode: str = 'periodization'):
434
"""
435
nD multiresolution analysis.
436
437
Parameters:
438
- data: Input nD array
439
- wavelet: Wavelet specification
440
- level: Number of decomposition levels (default: maximum possible)
441
- axes: Axes along which to perform transform (default: all axes)
442
- transform: Transform type ('swtn' for nD stationary WT)
443
- mode: Signal extension mode (default: 'periodization')
444
445
Returns:
446
List of nD approximation components at different scales
447
"""
448
449
def imran(mra_coeffs):
450
"""
451
Inverse nD MRA via summation.
452
453
Parameters:
454
- mra_coeffs: List of nD approximation components from mran()
455
456
Returns:
457
Reconstructed nD array (sum of all approximation levels)
458
"""
459
```
460
461
#### Usage Examples
462
463
```python
464
import pywt
465
import numpy as np
466
467
# 3D volume multiresolution analysis
468
volume_shape = (64, 64, 64)
469
x, y, z = np.mgrid[0:volume_shape[0], 0:volume_shape[1], 0:volume_shape[2]]
470
471
# Create 3D test volume with multiple scales
472
volume = (0.5 * np.sin(2 * np.pi * x / 32) * np.cos(2 * np.pi * y / 32) * np.sin(2 * np.pi * z / 32) +
473
0.3 * np.sin(2 * np.pi * x / 8) * np.cos(2 * np.pi * y / 8) * np.sin(2 * np.pi * z / 8) +
474
0.1 * np.random.randn(*volume_shape))
475
476
print(f"3D volume shape: {volume.shape}")
477
478
# 3D MRA
479
mra3d_coeffs = pywt.mran(volume, 'haar', level=4, transform='swtn')
480
print(f"3D MRA levels: {len(mra3d_coeffs)}")
481
print(f"3D approximation shapes: {[approx.shape for approx in mra3d_coeffs]}")
482
483
# Perfect reconstruction
484
reconstructed_3d = pywt.imran(mra3d_coeffs)
485
reconstruction_error_3d = np.max(np.abs(volume - reconstructed_3d))
486
print(f"3D MRA reconstruction error: {reconstruction_error_3d:.2e}")
487
488
# Analyze 3D scales
489
print("\n3D Scale Analysis:")
490
for i, component in enumerate(mra3d_coeffs):
491
scale = 2**(i+1)
492
energy = np.sum(component**2)
493
mean_val = np.mean(component)
494
std_val = np.std(component)
495
print(f"Level {i+1} (scale {scale}): energy = {energy:.2f}, mean = {mean_val:.4f}, std = {std_val:.4f}")
496
497
# Time series MRA (1D signal in higher dimensional context)
498
time_series_length = 10000
499
time_series = np.cumsum(np.random.randn(time_series_length)) * 0.1 # Random walk
500
time_series += np.sin(2 * np.pi * np.arange(time_series_length) / 100) # Add trend
501
502
# Reshape as pseudo-2D for demonstration
503
ts_reshaped = time_series.reshape(100, 100)
504
mra_ts = pywt.mra2(ts_reshaped, 'db8', level=6)
505
506
print(f"\nTime series MRA:")
507
print(f"Original time series length: {time_series_length}")
508
print(f"Reshaped as: {ts_reshaped.shape}")
509
print(f"MRA levels: {len(mra_ts)}")
510
511
# Reconstruct and verify
512
reconstructed_ts = pywt.imra2(mra_ts)
513
ts_error = np.max(np.abs(ts_reshaped - reconstructed_ts))
514
print(f"Time series MRA error: {ts_error:.2e}")
515
516
# Demonstrate memory efficiency of MRA vs full wavelet decomposition
517
def compare_memory_usage(data_shape, wavelet, level):
518
"""Compare memory usage between MRA and full wavelet decomposition."""
519
520
# Create dummy data
521
data = np.random.randn(*data_shape)
522
523
# MRA memory (approximations only)
524
mra_result = pywt.mra(data.ravel(), wavelet, level=level) if len(data_shape) == 1 else \
525
pywt.mra2(data, wavelet, level=level) if len(data_shape) == 2 else \
526
pywt.mran(data, wavelet, level=level)
527
528
mra_memory = sum(component.nbytes for component in mra_result)
529
530
# Full decomposition memory (approximation + all details)
531
if len(data_shape) == 1:
532
full_coeffs = pywt.wavedec(data.ravel(), wavelet, level=level)
533
full_memory = sum(coeff.nbytes for coeff in full_coeffs)
534
elif len(data_shape) == 2:
535
full_coeffs = pywt.wavedec2(data, wavelet, level=level)
536
full_memory = full_coeffs[0].nbytes + sum(sum(detail.nbytes for detail in level_details)
537
for level_details in full_coeffs[1:])
538
else:
539
full_coeffs = pywt.wavedecn(data, wavelet, level=level)
540
full_memory = full_coeffs[0].nbytes + sum(sum(detail.nbytes for detail in level_details.values())
541
for level_details in full_coeffs[1:])
542
543
original_memory = data.nbytes
544
545
return {
546
'original': original_memory,
547
'mra': mra_memory,
548
'full': full_memory,
549
'mra_ratio': mra_memory / original_memory,
550
'full_ratio': full_memory / original_memory
551
}
552
553
# Memory comparison for different data sizes
554
test_shapes = [(2048,), (512, 512), (64, 64, 64)]
555
556
print("\nMemory Usage Comparison:")
557
print(f"{'Shape':<15} {'Original (MB)':<12} {'MRA (MB)':<10} {'Full (MB)':<10} {'MRA Ratio':<10} {'Full Ratio':<10}")
558
print("-" * 80)
559
560
for shape in test_shapes:
561
memory_info = compare_memory_usage(shape, 'db4', 5)
562
print(f"{str(shape):<15} {memory_info['original']/1024**2:<12.2f} "
563
f"{memory_info['mra']/1024**2:<10.2f} {memory_info['full']/1024**2:<10.2f} "
564
f"{memory_info['mra_ratio']:<10.2f} {memory_info['full_ratio']:<10.2f}")
565
```
566
567
## Types
568
569
```python { .api }
570
# MRA transform backends
571
MRATransform1D = Literal['swt', 'dwt']
572
MRATransform2D = Literal['swt2', 'dwt2']
573
MRATransformND = Literal['swtn', 'dwtn']
574
575
# MRA coefficient format (list of approximation components)
576
MRACoeffs1D = List[np.ndarray] # [A1, A2, ..., An] - finest to coarsest
577
MRACoeffs2D = List[np.ndarray] # [A1, A2, ..., An] - 2D approximations
578
MRACoeffsND = List[np.ndarray] # [A1, A2, ..., An] - nD approximations
579
```