0
# Continuous Wavelet Transform
1

2
Continuous wavelet transform (CWT) for time-frequency analysis providing detailed spectral analysis with adjustable time-frequency resolution for non-stationary signal analysis.
3

4
## Capabilities
5

6
### Continuous Wavelet Transform
7

8
Time-frequency decomposition using continuous wavelets with arbitrary scales.
9

10
```python { .api }
11
def cwt(data, scales, wavelet, sampling_period: float = 1.0, 
12
        method: str = 'conv', axis: int = -1, *, precision: int = 12):
13
    """
14
    1D continuous wavelet transform.
15
    
16
    Parameters:
17
    - data: Input 1D signal
18
    - scales: Array of scales for the transform
19
    - wavelet: Continuous wavelet name or ContinuousWavelet object
20
    - sampling_period: Sampling period for frequency calculation (default: 1.0)
21
    - method: Computation method ('conv' for convolution, 'fft' for FFT-based)
22
    - axis: Axis along which to perform CWT (default: -1)
23
    - precision: Precision for wavelet function approximation (default: 12)
24
    
25
    Returns:
26
    (coefficients, frequencies) tuple where:
27
    - coefficients: 2D array (len(scales), len(data)) of CWT coefficients
28
    - frequencies: Array of frequencies corresponding to scales
29
    """
30
```
31

32
#### Usage Examples
33

34
```python
35
import pywt
36
import numpy as np
37
import matplotlib.pyplot as plt
38

39
# Create test signal with time-varying frequency
40
dt = 0.01
41
t = np.arange(0, 10, dt)
42
freq1, freq2 = 2, 10
43

44
# Chirp signal: frequency increases linearly with time
45
signal = np.sin(2 * np.pi * (freq1 + (freq2 - freq1) * t / max(t)) * t)
46

47
# Add a transient high-frequency component
48
transient_start, transient_end = 300, 400  # Sample indices
49
signal[transient_start:transient_end] += 2 * np.sin(2 * np.pi * 50 * t[transient_start:transient_end])
50

51
# Add noise
52
noise_level = 0.3
53
noisy_signal = signal + noise_level * np.random.randn(len(signal))
54

55
print(f"Signal length: {len(noisy_signal)}")
56
print(f"Sampling period: {dt}")
57
print(f"Total duration: {t[-1]:.2f} seconds")
58

59
# Define scales for CWT
60
scales = np.arange(1, 128)  # Logarithmic scale often better
61
# scales = np.logspace(0, 2, 64)  # Alternative: logarithmic spacing
62

63
# Perform CWT using different wavelets
64
wavelets = ['mexh', 'morl', 'cgau8']
65
cwt_results = {}
66

67
for wavelet_name in wavelets:
68
    coefficients, frequencies = pywt.cwt(noisy_signal, scales, wavelet_name, sampling_period=dt)
69
    cwt_results[wavelet_name] = (coefficients, frequencies)
70
    print(f"{wavelet_name}: CWT shape {coefficients.shape}, frequency range {frequencies.min():.3f}-{frequencies.max():.3f} Hz")
71

72
# Visualization
73
fig, axes = plt.subplots(len(wavelets) + 1, 1, figsize=(15, 12))
74

75
# Original signal
76
axes[0].plot(t, signal, 'b-', label='Clean signal', alpha=0.7)
77
axes[0].plot(t, noisy_signal, 'r-', label='Noisy signal', alpha=0.8)
78
axes[0].set_title('Original Signal')
79
axes[0].set_xlabel('Time (s)')
80
axes[0].set_ylabel('Amplitude')
81
axes[0].legend()
82
axes[0].grid(True)
83

84
# CWT plots for each wavelet
85
for i, wavelet_name in enumerate(wavelets):
86
    coefficients, frequencies = cwt_results[wavelet_name]
87
    
88
    # Create time-frequency plot
89
    T, F = np.meshgrid(t, frequencies)
90
    im = axes[i+1].contourf(T, F, np.abs(coefficients), levels=50, cmap='jet')
91
    axes[i+1].set_title(f'CWT Scalogram - {wavelet_name.upper()}')
92
    axes[i+1].set_xlabel('Time (s)')
93
    axes[i+1].set_ylabel('Frequency (Hz)')
94
    plt.colorbar(im, ax=axes[i+1], label='|CWT Coefficients|')
95

96
plt.tight_layout()
97
plt.show()
98

99
# Advanced analysis: Ridge detection for instantaneous frequency
100
def detect_ridges(coefficients, frequencies, threshold_factor=0.5):
101
    """Detect ridges in CWT scalogram for instantaneous frequency estimation."""
102
    abs_coeffs = np.abs(coefficients)
103
    threshold = threshold_factor * np.max(abs_coeffs)
104
    
105
    ridges = []
106
    for t_idx in range(abs_coeffs.shape[1]):
107
        column = abs_coeffs[:, t_idx]
108
        peaks = []
109
        for f_idx in range(1, len(column)-1):
110
            if (column[f_idx] > column[f_idx-1] and 
111
                column[f_idx] > column[f_idx+1] and
112
                column[f_idx] > threshold):
113
                peaks.append((f_idx, column[f_idx]))
114
        
115
        if peaks:
116
            # Take the strongest peak
117
            strongest_peak = max(peaks, key=lambda x: x[1])
118
            ridges.append(frequencies[strongest_peak[0]])
119
        else:
120
            ridges.append(np.nan)
121
    
122
    return np.array(ridges)
123

124
# Detect instantaneous frequency using Morlet wavelet
125
coefficients, frequencies = cwt_results['morl']
126
inst_freq = detect_ridges(coefficients, frequencies)
127

128
# Plot instantaneous frequency
129
plt.figure(figsize=(12, 6))
130
plt.subplot(2, 1, 1)
131
plt.plot(t, noisy_signal)
132
plt.title('Signal')
133
plt.ylabel('Amplitude')
134

135
plt.subplot(2, 1, 2)
136
plt.plot(t, inst_freq, 'r-', linewidth=2, label='Estimated Instantaneous Frequency')
137
# True instantaneous frequency for the chirp
138
true_inst_freq = freq1 + (freq2 - freq1) * t / max(t)
139
plt.plot(t, true_inst_freq, 'b--', linewidth=2, label='True Instantaneous Frequency')
140
plt.xlabel('Time (s)')
141
plt.ylabel('Frequency (Hz)')
142
plt.legend()
143
plt.title('Instantaneous Frequency Estimation')
144
plt.grid(True)
145
plt.tight_layout()
146
plt.show()
147

148
# CWT-based denoising using thresholding
149
def cwt_denoise(signal, scales, wavelet, threshold_factor=0.1):
150
    """Denoise signal using CWT thresholding."""
151
    coefficients, frequencies = pywt.cwt(signal, scales, wavelet)
152
    
153
    # Threshold coefficients
154
    threshold = threshold_factor * np.max(np.abs(coefficients))
155
    coefficients_thresh = np.where(np.abs(coefficients) > threshold, coefficients, 0)
156
    
157
    # Inverse CWT (approximate reconstruction using mother wavelet)
158
    # Note: Exact inverse CWT requires admissibility constant
159
    reconstruction = np.zeros_like(signal)
160
    
161
    wavelet_obj = pywt.ContinuousWavelet(wavelet)
162
    _, psi = wavelet_obj.wavefun()
163
    
164
    for i, scale in enumerate(scales):
165
        # Simple reconstruction (not exact inverse CWT)
166
        contribution = np.real(coefficients_thresh[i, :])
167
        reconstruction += contribution / scale
168
    
169
    return reconstruction / len(scales)
170

171
# Apply CWT denoising
172
denoised_cwt = cwt_denoise(noisy_signal, scales[:32], 'morl', threshold_factor=0.2)
173

174
plt.figure(figsize=(12, 6))
175
plt.plot(t, signal, 'b-', label='Original clean signal', linewidth=2)
176
plt.plot(t, noisy_signal, 'r-', alpha=0.7, label='Noisy signal')
177
plt.plot(t, denoised_cwt, 'g-', label='CWT denoised', linewidth=2)
178
plt.xlabel('Time (s)')
179
plt.ylabel('Amplitude')
180
plt.title('CWT-based Denoising')
181
plt.legend()
182
plt.grid(True)
183
plt.show()
184

185
# SNR calculation
186
def calculate_snr(clean, noisy):
187
    """Calculate Signal-to-Noise Ratio in dB."""
188
    signal_power = np.mean(clean**2)
189
    noise_power = np.mean((noisy - clean)**2)
190
    return 10 * np.log10(signal_power / noise_power)
191

192
original_snr = calculate_snr(signal, noisy_signal)
193
denoised_snr = calculate_snr(signal, denoised_cwt)
194
print(f"Original SNR: {original_snr:.2f} dB")
195
print(f"Denoised SNR: {denoised_snr:.2f} dB")
196
print(f"SNR improvement: {denoised_snr - original_snr:.2f} dB")
197
```
198

199
### Scale-Frequency Conversion Utilities
200

201
Utility functions for converting between CWT scales and frequencies.
202

203
```python { .api }
204
def scale2frequency(wavelet, scale: float, precision: int = 8) -> float:
205
    """
206
    Convert CWT scale to normalized frequency.
207
    
208
    Parameters:
209
    - wavelet: Continuous wavelet specification
210
    - scale: CWT scale value
211
    - precision: Precision for wavelet function approximation
212
    
213
    Returns:
214
    Normalized frequency (sampling frequency = 1.0)
215
    """
216

217
def frequency2scale(wavelet, freq: float, precision: int = 8) -> float:
218
    """
219
    Convert normalized frequency to CWT scale.
220
    
221
    Parameters:
222
    - wavelet: Continuous wavelet specification  
223
    - freq: Normalized frequency
224
    - precision: Precision for wavelet function approximation
225
    
226
    Returns:
227
    CWT scale value
228
    """
229

230
def next_fast_len(n: int) -> int:
231
    """
232
    Round up size to the nearest power of two for FFT optimization.
233
    
234
    Parameters:
235
    - n: Input size
236
    
237
    Returns:
238
    Next power of two >= n for efficient FFT computation
239
    """
240
```
241

242
#### Usage Examples
243

244
```python
245
import pywt
246
import numpy as np
247
import matplotlib.pyplot as plt
248

249
# Scale-frequency relationship analysis
250
wavelets = ['mexh', 'morl', 'cgau8', 'cmor1.5-1.0']
251
scales = np.arange(1, 101)
252

253
plt.figure(figsize=(12, 8))
254

255
for wavelet_name in wavelets:
256
    try:
257
        frequencies = []
258
        for scale in scales:
259
            freq = pywt.scale2frequency(wavelet_name, scale)
260
            frequencies.append(freq)
261
        
262
        plt.loglog(scales, frequencies, 'o-', label=wavelet_name, markersize=3)
263
        
264
        # Show inverse relationship  
265
        print(f"{wavelet_name}: Scale 10 -> Frequency {pywt.scale2frequency(wavelet_name, 10):.4f}")
266
        freq_01 = 0.1
267
        scale_back = pywt.frequency2scale(wavelet_name, freq_01)
268
        print(f"{wavelet_name}: Frequency 0.1 -> Scale {scale_back:.4f}")
269
        
270
    except Exception as e:
271
        print(f"Error with {wavelet_name}: {e}")
272

273
plt.xlabel('Scale')
274
plt.ylabel('Normalized Frequency')
275
plt.title('Scale-Frequency Relationship for Different Wavelets')
276
plt.legend()
277
plt.grid(True, which="both", ls="-", alpha=0.3)
278
plt.show()
279

280
# Practical example: Design CWT analysis for specific frequency range
281
target_freq_range = [1, 50]  # Hz
282
sampling_freq = 1000  # Hz
283
sampling_period = 1.0 / sampling_freq
284

285
# Convert to normalized frequencies
286
norm_freq_range = [f * sampling_period for f in target_freq_range]
287
print(f"Target frequency range: {target_freq_range} Hz")
288
print(f"Normalized frequency range: {norm_freq_range}")
289

290
# Find corresponding scales for Morlet wavelet
291
wavelet = 'morl'
292
scales_for_range = []
293
for norm_freq in norm_freq_range:
294
    scale = pywt.frequency2scale(wavelet, norm_freq)
295
    scales_for_range.append(scale)
296

297
print(f"Corresponding scales for {wavelet}: {scales_for_range}")
298

299
# Create logarithmically spaced scales covering the range
300
num_scales = 64
301
scale_min, scale_max = min(scales_for_range), max(scales_for_range)
302
scales_log = np.logspace(np.log10(scale_min), np.log10(scale_max), num_scales)
303

304
# Verify frequency coverage
305
freq_coverage = [pywt.scale2frequency(wavelet, s) / sampling_period for s in scales_log]
306
print(f"Frequency coverage: {min(freq_coverage):.2f} - {max(freq_coverage):.2f} Hz")
307

308
# Test with synthetic signal in target range
309
t = np.arange(0, 2, sampling_period)
310
test_signal = (np.sin(2 * np.pi * 5 * t) +  # 5 Hz
311
               np.sin(2 * np.pi * 25 * t))   # 25 Hz
312

313
coefficients, frequencies = pywt.cwt(test_signal, scales_log, wavelet, sampling_period=sampling_period)
314

315
# Plot result
316
plt.figure(figsize=(12, 6))
317
T, F = np.meshgrid(t, frequencies)
318
plt.contourf(T, F, np.abs(coefficients), levels=50, cmap='jet')
319
plt.colorbar(label='|CWT Coefficients|')
320
plt.xlabel('Time (s)')
321
plt.ylabel('Frequency (Hz)')
322
plt.title('CWT Analysis in Target Frequency Range')
323
plt.ylim(target_freq_range)
324
plt.show()
325
```
326

327
### Complex CWT Analysis
328

329
Working with complex-valued continuous wavelets for phase and amplitude analysis.
330

331
#### Usage Examples
332

333
```python
334
import pywt
335
import numpy as np
336
import matplotlib.pyplot as plt
337

338
# Create test signal with amplitude and frequency modulation
339
t = np.linspace(0, 4, 1000)
340
dt = t[1] - t[0]
341

342
# Amplitude modulated signal
343
carrier_freq = 10  # Hz
344
mod_freq = 1       # Hz
345
am_signal = (1 + 0.5 * np.sin(2 * np.pi * mod_freq * t)) * np.sin(2 * np.pi * carrier_freq * t)
346

347
# Frequency modulated signal  
348
freq_deviation = 5  # Hz
349
fm_signal = np.sin(2 * np.pi * (carrier_freq + freq_deviation * np.sin(2 * np.pi * mod_freq * t)) * t)
350

351
# Combine signals in different time windows
352
signal = np.concatenate([am_signal[:250], fm_signal[250:500], 
353
                        am_signal[500:750], fm_signal[750:]])
354

355
print(f"Signal length: {len(signal)}")
356

357
# CWT with complex Morlet wavelet
358
scales = np.arange(1, 64)
359
coefficients, frequencies = pywt.cwt(signal, scales, 'morl', sampling_period=dt)
360

361
print(f"CWT coefficients are complex: {np.iscomplexobj(coefficients)}")
362
print(f"CWT shape: {coefficients.shape}")
363

364
# Extract amplitude and phase
365
amplitude = np.abs(coefficients)
366
phase = np.angle(coefficients)
367
instantaneous_phase = np.unwrap(phase, axis=1)
368

369
# Plot complex CWT analysis
370
fig, axes = plt.subplots(4, 1, figsize=(15, 12))
371

372
# Original signal
373
axes[0].plot(t, signal)
374
axes[0].set_title('Test Signal (AM + FM)')
375
axes[0].set_ylabel('Amplitude')
376
axes[0].grid(True)
377

378
# Amplitude scalogram
379
T, F = np.meshgrid(t, frequencies)
380
im1 = axes[1].contourf(T, F, amplitude, levels=50, cmap='jet')
381
axes[1].set_title('CWT Amplitude Scalogram')
382
axes[1].set_ylabel('Frequency (Hz)')
383
plt.colorbar(im1, ax=axes[1], label='Amplitude')
384

385
# Phase scalogram
386
im2 = axes[2].contourf(T, F, phase, levels=50, cmap='hsv')
387
axes[2].set_title('CWT Phase Scalogram') 
388
axes[2].set_ylabel('Frequency (Hz)')
389
plt.colorbar(im2, ax=axes[2], label='Phase (rad)')
390

391
# Instantaneous frequency (derivative of phase)
392
# Focus on carrier frequency region
393
carrier_idx = np.argmin(np.abs(frequencies - carrier_freq))
394
inst_freq_approx = np.diff(instantaneous_phase[carrier_idx, :]) / (2 * np.pi * dt)
395

396
axes[3].plot(t[:-1], inst_freq_approx, 'r-', linewidth=2, label='Estimated Inst. Freq.')
397
axes[3].axhline(y=carrier_freq, color='b', linestyle='--', label=f'Carrier ({carrier_freq} Hz)')
398
axes[3].set_title('Instantaneous Frequency Estimation')
399
axes[3].set_xlabel('Time (s)')
400
axes[3].set_ylabel('Frequency (Hz)')
401
axes[3].legend()
402
axes[3].grid(True)
403

404
plt.tight_layout()
405
plt.show()
406

407
# Amplitude and phase tracking at specific frequency
408
target_freq = carrier_freq
409
freq_idx = np.argmin(np.abs(frequencies - target_freq))
410

411
amplitude_track = amplitude[freq_idx, :]
412
phase_track = phase[freq_idx, :]
413

414
plt.figure(figsize=(12, 8))
415

416
plt.subplot(3, 1, 1)
417
plt.plot(t, signal)
418
plt.title(f'Original Signal')
419
plt.ylabel('Amplitude')
420
plt.grid(True)
421

422
plt.subplot(3, 1, 2)
423
plt.plot(t, amplitude_track, 'r-', linewidth=2)
424
plt.title(f'Amplitude Tracking at {target_freq} Hz')
425
plt.ylabel('CWT Amplitude')
426
plt.grid(True)
427

428
plt.subplot(3, 1, 3)
429
plt.plot(t, phase_track, 'g-', linewidth=2)
430
plt.title(f'Phase Tracking at {target_freq} Hz')
431
plt.xlabel('Time (s)')
432
plt.ylabel('Phase (rad)')
433
plt.grid(True)
434

435
plt.tight_layout()
436
plt.show()
437

438
# Demonstrate admissibility of different wavelets
439
wavelets_to_test = ['mexh', 'morl', 'cgau8', 'cmor1.5-1.0']
440

441
for wavelet_name in wavelets_to_test:
442
    try:
443
        wavelet_obj = pywt.ContinuousWavelet(wavelet_name)
444
        print(f"\n{wavelet_name.upper()}:")
445
        print(f"  Complex CWT: {wavelet_obj.complex_cwt}")
446
        print(f"  Center frequency: {wavelet_obj.center_frequency:.4f}")
447
        print(f"  Bandwidth frequency: {wavelet_obj.bandwidth_frequency:.4f}")
448
        print(f"  Support: [{wavelet_obj.lower_bound:.2f}, {wavelet_obj.upper_bound:.2f}]")
449
        
450
        # Show wavelet function
451
        psi, x = wavelet_obj.wavefun()
452
        
453
        if wavelet_obj.complex_cwt:
454
            print(f"  Real part range: [{np.min(np.real(psi)):.4f}, {np.max(np.real(psi)):.4f}]")
455
            print(f"  Imag part range: [{np.min(np.imag(psi)):.4f}, {np.max(np.imag(psi)):.4f}]")
456
        else:
457
            print(f"  Value range: [{np.min(psi):.4f}, {np.max(psi):.4f}]")
458
            
459
    except Exception as e:
460
        print(f"Error with {wavelet_name}: {e}")
461
```
462

463
## Types
464

465
```python { .api }
466
# CWT result format
467
CWTResult = Tuple[np.ndarray, np.ndarray]  # (coefficients, frequencies)
468

469
# CWT coefficients are 2D: (len(scales), len(data))
470
# Can be real or complex depending on wavelet
471
CWTCoefficients = Union[np.ndarray[np.float64], np.ndarray[np.complex128]]
472

473
# Scales specification
474
Scales = Union[np.ndarray, list, range]
475

476
# CWT methods
477
CWTMethod = Literal['conv', 'fft']
478

479
# Common continuous wavelets
480
ContinuousWaveletName = Literal[
481
    'mexh',           # Mexican hat
482
    'morl',           # Morlet 
483
    'cgau1', 'cgau2', 'cgau3', 'cgau4', 'cgau5', 'cgau6', 'cgau7', 'cgau8',  # Complex Gaussian
484
    'shan',           # Shannon
485
    'fbsp',           # Frequency B-Spline
486
    'cmor',           # Complex Morlet (cmor<bandwidth>-<center_frequency>)
487
    'gaus1', 'gaus2', 'gaus3', 'gaus4', 'gaus5', 'gaus6', 'gaus7', 'gaus8'   # Gaussian derivatives
488
]
489
```