0
# Correlation and Signal Processing
1

2
Correlation analysis, signal generation, and essential signal processing utilities for spectral analysis applications. These functions provide the fundamental building blocks for many spectral estimation methods and signal analysis tasks.
3

4
## Capabilities
5

6
### Correlation Functions
7

8
Computation of auto-correlation and cross-correlation functions with various normalization options.
9

10
```python { .api }
11
def CORRELATION(x, y=None, maxlags=None, norm='unbiased'):
12
    """
13
    Correlation function computation (auto or cross-correlation).
14
    
15
    Parameters:
16
    - x: array-like, first input signal
17
    - y: array-like, second input signal (None for autocorrelation)
18
    - maxlags: int, maximum lag to compute (None for full range)
19
    - norm: str, normalization type ('biased', 'unbiased', 'coeff', None)
20
    
21
    Returns:
22
    array: Correlation values for positive lags only
23
    """
24

25
def xcorr(x, y=None, maxlags=None, norm='biased'):
26
    """
27
    Fast cross-correlation using numpy.correlate.
28
    
29
    Parameters:
30
    - x: array-like, first input signal
31
    - y: array-like, second input signal (None for autocorrelation)
32
    - maxlags: int, maximum lag to compute
33
    - norm: str, normalization type ('biased', 'unbiased', 'coeff', None)
34
    
35
    Returns:
36
    tuple: (correlation sequence array, lags array)
37
    """
38
```
39

40
### Signal Generation
41

42
Functions for generating test signals, waveforms, and synthetic data for spectral analysis.
43

44
```python { .api }
45
def chirp(t, f0=0., t1=1., f1=100., form='linear', phase=0):
46
    """
47
    Generate chirp (frequency-swept) signal.
48
    
49
    Parameters:
50
    - t: array-like, time vector
51
    - f0: float, starting frequency
52
    - t1: float, end time for frequency sweep
53
    - f1: float, ending frequency
54
    - form: str, sweep type ('linear', 'quadratic', 'logarithmic', 'hyperbolic')
55
    - phase: float, initial phase in radians
56
    
57
    Returns:
58
    array: Generated chirp signal
59
    """
60

61
def morlet(lb, ub, n):
62
    """
63
    Generate Morlet wavelet for time-frequency analysis.
64
    
65
    Parameters:
66
    - lb: float, lower bound of time interval
67
    - ub: float, upper bound of time interval
68
    - n: int, number of samples
69
    
70
    Returns:
71
    array: Morlet wavelet samples
72
    """
73

74
def mexican(lb, ub, n):
75
    """
76
    Generate Mexican hat wavelet (second derivative of Gaussian).
77
    
78
    Parameters:
79
    - lb: float, lower bound of time interval
80
    - ub: float, upper bound of time interval  
81
    - n: int, number of samples
82
    
83
    Returns:
84
    array: Mexican hat wavelet samples
85
    """
86

87
def meyeraux(x):
88
    """
89
    Meyer wavelet auxiliary function.
90
    
91
    Parameters:
92
    - x: array-like, input values
93
    
94
    Returns:
95
    array: Meyer auxiliary function values
96
    """
97
```
98

99
### Power and Decibel Conversions
100

101
Essential conversion functions between linear power, magnitude, and decibel scales.
102

103
```python { .api }
104
def pow2db(x):
105
    """
106
    Convert power values to decibels.
107
    
108
    Parameters:
109
    - x: array-like, power values (linear scale)
110
    
111
    Returns:
112
    array: Power in decibels (10*log10(x))
113
    """
114

115
def db2pow(xdb):
116
    """
117
    Convert decibel values to linear power.
118
    
119
    Parameters:
120
    - xdb: array-like, power values in dB
121
    
122
    Returns:
123
    array: Linear power values (10^(xdb/10))
124
    """
125

126
def mag2db(x):
127
    """
128
    Convert magnitude values to decibels.
129
    
130
    Parameters:
131
    - x: array-like, magnitude values
132
    
133
    Returns:
134
    array: Magnitude in decibels (20*log10(x))
135
    """
136

137
def db2mag(xdb):
138
    """
139
    Convert decibel values to linear magnitude.
140
    
141
    Parameters:
142
    - xdb: array-like, magnitude values in dB
143
    
144
    Returns:
145
    array: Linear magnitude values (10^(xdb/20))
146
    """
147
```
148

149
### Utility Functions
150

151
Mathematical and signal processing utility functions for various applications.
152

153
```python { .api }
154
def nextpow2(x):
155
    """
156
    Find the next power of 2 greater than or equal to x.
157
    
158
    Parameters:
159
    - x: float or int, input value
160
    
161
    Returns:
162
    int: Next power of 2
163
    """
164

165
def fftshift(x):
166
    """
167
    Wrapper to numpy.fft.fftshift for frequency domain centering.
168
    
169
    Parameters:
170
    - x: array-like, input array
171
    
172
    Returns:
173
    array: Shifted array with zero frequency at center
174
    """
175

176
def cshift(data, offset):
177
    """
178
    Circular shift of array elements.
179
    
180
    Parameters:
181
    - data: array-like, input data array
182
    - offset: int, number of positions to shift
183
    
184
    Returns:
185
    array: Circularly shifted array
186
    """
187
```
188

189
### Data Format Conversions
190

191
Functions for converting between different PSD representation formats.
192

193
```python { .api }
194
def twosided_2_onesided(data):
195
    """
196
    Convert two-sided PSD to one-sided PSD.
197
    
198
    Parameters:
199
    - data: array-like, two-sided PSD values
200
    
201
    Returns:
202
    array: One-sided PSD values (positive frequencies only)
203
    """
204

205
def onesided_2_twosided(data):
206
    """
207
    Convert one-sided PSD to two-sided PSD.
208
    
209
    Parameters:
210
    - data: array-like, one-sided PSD values
211
    
212
    Returns:
213
    array: Two-sided PSD values (negative and positive frequencies)
214
    """
215

216
def twosided_2_centerdc(data):
217
    """
218
    Convert two-sided PSD to center-DC format.
219
    
220
    Parameters:
221
    - data: array-like, two-sided PSD values
222
    
223
    Returns:
224
    array: Center-DC format PSD (DC at center)
225
    """
226

227
def centerdc_2_twosided(data):
228
    """
229
    Convert center-DC format PSD to standard two-sided format.
230
    
231
    Parameters:
232
    - data: array-like, center-DC format PSD
233
    
234
    Returns:
235
    array: Two-sided PSD values
236
    """
237
```
238

239
## Usage Examples
240

241
### Correlation Analysis
242

243
```python
244
import spectrum
245
import numpy as np
246
import matplotlib.pyplot as plt
247

248
# Generate test signals for correlation analysis
249
N = 200
250
n = np.arange(N)
251

252
# Signal 1: Sinusoid + noise
253
f1 = 0.1
254
signal1 = np.sin(2*np.pi*f1*n) + 0.3*np.random.randn(N)
255

256
# Signal 2: Delayed and noisy version of signal1
257
delay = 20
258
signal2 = np.zeros(N)
259
signal2[delay:] = signal1[:-delay] + 0.2*np.random.randn(N-delay)
260

261
plt.figure(figsize=(15, 12))
262

263
# Plot signals
264
plt.subplot(3, 2, 1)
265
plt.plot(signal1, 'b-', label='Signal 1', alpha=0.8)
266
plt.plot(signal2, 'r-', label='Signal 2 (delayed)', alpha=0.8)
267
plt.title('Input Signals')
268
plt.xlabel('Sample')
269
plt.ylabel('Amplitude')
270
plt.legend()
271
plt.grid(True)
272

273
# Autocorrelation of signal 1
274
plt.subplot(3, 2, 2)
275
maxlags = 50
276
autocorr1 = spectrum.CORRELATION(signal1, maxlags=maxlags, norm='unbiased')
277
lags_pos = np.arange(maxlags+1)
278
plt.plot(lags_pos, autocorr1, 'b-', linewidth=2)
279
plt.title('Autocorrelation of Signal 1')
280
plt.xlabel('Lag')
281
plt.ylabel('Correlation')
282
plt.grid(True)
283

284
# Cross-correlation using CORRELATION function
285
plt.subplot(3, 2, 3)
286
crosscorr = spectrum.CORRELATION(signal1, signal2, maxlags=maxlags, norm='unbiased')
287
plt.plot(lags_pos, crosscorr, 'g-', linewidth=2)
288
plt.title('Cross-correlation (CORRELATION function)')
289
plt.xlabel('Lag')
290
plt.ylabel('Correlation')
291
plt.grid(True)
292

293
# Cross-correlation using xcorr function (returns both positive and negative lags)
294
plt.subplot(3, 2, 4)
295
crosscorr_full, lags_full = spectrum.xcorr(signal1, signal2, maxlags=maxlags, norm='unbiased')
296
plt.plot(lags_full, crosscorr_full, 'r-', linewidth=2)
297
plt.axvline(delay, color='k', linestyle='--', alpha=0.7, label=f'True delay={delay}')
298
peak_lag = lags_full[np.argmax(crosscorr_full)]
299
plt.axvline(peak_lag, color='g', linestyle='--', alpha=0.7, label=f'Detected delay={peak_lag}')
300
plt.title('Cross-correlation (xcorr function)')
301
plt.xlabel('Lag')
302
plt.ylabel('Correlation')
303
plt.legend()
304
plt.grid(True)
305

306
# Compare different normalization methods
307
plt.subplot(3, 2, 5)
308
norms = ['biased', 'unbiased', 'coeff', None]
309
for norm_type in norms:
310
    if norm_type is not None:
311
        corr_norm = spectrum.CORRELATION(signal1, maxlags=30, norm=norm_type)
312
        label = f'{norm_type}'
313
    else:
314
        corr_norm = spectrum.CORRELATION(signal1, maxlags=30, norm=None)
315
        label = 'unnormalized'
316
    
317
    lags_norm = np.arange(len(corr_norm))
318
    plt.plot(lags_norm, corr_norm, linewidth=2, label=label)
319

320
plt.title('Autocorrelation: Different Normalizations')
321
plt.xlabel('Lag')
322
plt.ylabel('Correlation')
323
plt.legend()
324
plt.grid(True)
325

326
# Practical application: Time delay estimation
327
plt.subplot(3, 2, 6)
328
delays_to_test = range(5, 35, 5)
329
detected_delays = []
330
correlation_peaks = []
331

332
for test_delay in delays_to_test:
333
    # Create delayed version
334
    sig2_test = np.zeros(N)
335
    sig2_test[test_delay:] = signal1[:-test_delay] + 0.1*np.random.randn(N-test_delay)
336
    
337
    # Find peak in cross-correlation
338
    cc, lags = spectrum.xcorr(signal1, sig2_test, maxlags=maxlags)
339
    peak_idx = np.argmax(cc)
340
    detected_delay = lags[peak_idx]
341
    peak_value = cc[peak_idx]
342
    
343
    detected_delays.append(detected_delay)
344
    correlation_peaks.append(peak_value)
345

346
plt.plot(delays_to_test, detected_delays, 'bo-', linewidth=2, label='Detected')
347
plt.plot(delays_to_test, delays_to_test, 'r--', linewidth=2, label='True delay')
348
plt.title('Time Delay Estimation Accuracy')
349
plt.xlabel('True Delay (samples)')
350
plt.ylabel('Detected Delay (samples)')
351
plt.legend()
352
plt.grid(True)
353

354
plt.tight_layout()
355
plt.show()
356

357
# Print delay estimation results
358
print("Delay Estimation Results:")
359
print("-" * 30)
360
for true_del, det_del, peak in zip(delays_to_test, detected_delays, correlation_peaks):
361
    error = abs(true_del - det_del)
362
    print(f"True: {true_del:2d}, Detected: {det_del:2d}, Error: {error:2d}, Peak: {peak:.3f}")
363
```
364

365
### Signal Generation and Analysis
366

367
```python
368
import spectrum
369
import numpy as np
370
import matplotlib.pyplot as plt
371

372
# Time parameters
373
fs = 1000  # Sampling frequency
374
T = 2.0    # Duration
375
t = np.linspace(0, T, int(fs*T), False)
376

377
# Generate different types of signals
378
plt.figure(figsize=(15, 12))
379

380
# 1. Linear chirp
381
plt.subplot(3, 3, 1)
382
chirp_linear = spectrum.chirp(t, f0=50, t1=T, f1=200, form='linear')
383
plt.plot(t, chirp_linear)
384
plt.title('Linear Chirp (50-200 Hz)')
385
plt.xlabel('Time (s)')
386
plt.ylabel('Amplitude')
387
plt.grid(True)
388

389
# 2. Logarithmic chirp
390
plt.subplot(3, 3, 2)
391
chirp_log = spectrum.chirp(t, f0=50, t1=T, f1=200, form='logarithmic')
392
plt.plot(t, chirp_log)
393
plt.title('Logarithmic Chirp (50-200 Hz)')
394
plt.xlabel('Time (s)')
395
plt.ylabel('Amplitude')
396
plt.grid(True)
397

398
# 3. Quadratic chirp
399
plt.subplot(3, 3, 3)
400
chirp_quad = spectrum.chirp(t, f0=50, t1=T, f1=200, form='quadratic')
401
plt.plot(t, chirp_quad)
402
plt.title('Quadratic Chirp (50-200 Hz)')
403
plt.xlabel('Time (s)')
404
plt.ylabel('Amplitude')
405
plt.grid(True)
406

407
# 4. Morlet wavelet
408
plt.subplot(3, 3, 4)
409
morlet_wave = spectrum.morlet(-3, 3, 200)
410
t_morlet = np.linspace(-3, 3, 200)
411
plt.plot(t_morlet, morlet_wave, 'r-', linewidth=2)
412
plt.title('Morlet Wavelet')
413
plt.xlabel('Time')
414
plt.ylabel('Amplitude')
415
plt.grid(True)
416

417
# 5. Mexican hat wavelet
418
plt.subplot(3, 3, 5)
419
mexican_wave = spectrum.mexican(-3, 3, 200)
420
plt.plot(t_morlet, mexican_wave, 'g-', linewidth=2)
421
plt.title('Mexican Hat Wavelet')
422
plt.xlabel('Time')
423
plt.ylabel('Amplitude')
424
plt.grid(True)
425

426
# Spectrograms of chirp signals
427
plt.subplot(3, 3, 7)
428
from scipy import signal
429
f_spec, t_spec, Sxx = signal.spectrogram(chirp_linear, fs, nperseg=256, noverlap=128)
430
plt.pcolormesh(t_spec, f_spec, 10*np.log10(Sxx), shading='gouraud')
431
plt.ylabel('Frequency (Hz)')
432
plt.xlabel('Time (s)')
433
plt.title('Linear Chirp Spectrogram')
434
plt.ylim(0, 300)
435

436
plt.subplot(3, 3, 8)
437
f_spec, t_spec, Sxx = signal.spectrogram(chirp_log, fs, nperseg=256, noverlap=128)
438
plt.pcolormesh(t_spec, f_spec, 10*np.log10(Sxx), shading='gouraud')
439
plt.ylabel('Frequency (Hz)')
440
plt.xlabel('Time (s)')
441
plt.title('Logarithmic Chirp Spectrogram')
442
plt.ylim(0, 300)
443

444
plt.subplot(3, 3, 9)
445
f_spec, t_spec, Sxx = signal.spectrogram(chirp_quad, fs, nperseg=256, noverlap=128)
446
plt.pcolormesh(t_spec, f_spec, 10*np.log10(Sxx), shading='gouraud')
447
plt.ylabel('Frequency (Hz)')
448
plt.xlabel('Time (s)')
449
plt.title('Quadratic Chirp Spectrogram')
450
plt.ylim(0, 300)
451

452
plt.tight_layout()
453
plt.show()
454

455
# Compare chirp types in frequency domain
456
plt.figure(figsize=(12, 8))
457

458
signals = {
459
    'Linear': chirp_linear,
460
    'Logarithmic': chirp_log, 
461
    'Quadratic': chirp_quad
462
}
463

464
for i, (chirp_type, chirp_signal) in enumerate(signals.items()):
465
    plt.subplot(2, 2, i+1)
466
    
467
    # Compute PSD using Welch method
468
    f_psd, psd = signal.welch(chirp_signal, fs, nperseg=1024, noverlap=512)
469
    plt.semilogy(f_psd, psd, linewidth=2)
470
    plt.title(f'{chirp_type} Chirp - PSD')
471
    plt.xlabel('Frequency (Hz)')
472
    plt.ylabel('PSD')
473
    plt.grid(True)
474
    plt.xlim(0, 300)
475

476
# Wavelet frequency content
477
plt.subplot(2, 2, 4)
478
morlet_fft = np.fft.fft(morlet_wave, 1024)
479
freqs_wav = np.fft.fftfreq(1024, d=6/200)  # Adjust for time span
480
plt.semilogy(np.fft.fftshift(freqs_wav), np.fft.fftshift(np.abs(morlet_fft)**2), 'r-', linewidth=2, label='Morlet')
481

482
mexican_fft = np.fft.fft(mexican_wave, 1024)
483
plt.semilogy(np.fft.fftshift(freqs_wav), np.fft.fftshift(np.abs(mexican_fft)**2), 'g-', linewidth=2, label='Mexican Hat')
484
plt.title('Wavelet Frequency Content')
485
plt.xlabel('Frequency')
486
plt.ylabel('Power')
487
plt.legend()
488
plt.grid(True)
489
plt.xlim(-2, 2)
490

491
plt.tight_layout()
492
plt.show()
493
```
494

495
### Power and Decibel Conversions
496

497
```python
498
import spectrum
499
import numpy as np
500
import matplotlib.pyplot as plt
501

502
# Generate test data with different dynamic ranges
503
linear_powers = np.logspace(-3, 3, 100)  # 10^-3 to 10^3
504
linear_magnitudes = np.sqrt(linear_powers)
505

506
# Convert to dB scales
507
power_db = spectrum.pow2db(linear_powers)
508
magnitude_db = spectrum.mag2db(linear_magnitudes)
509

510
# Convert back to linear
511
power_recovered = spectrum.db2pow(power_db)
512
magnitude_recovered = spectrum.db2mag(magnitude_db)
513

514
plt.figure(figsize=(15, 10))
515

516
# Linear vs dB power scales
517
plt.subplot(2, 3, 1)
518
plt.loglog(linear_powers, linear_powers, 'b-', linewidth=2, label='Linear scale')
519
plt.xlabel('Input Power')
520
plt.ylabel('Power')
521
plt.title('Linear Power Scale')
522
plt.grid(True)
523
plt.legend()
524

525
plt.subplot(2, 3, 2)
526
plt.semilogx(linear_powers, power_db, 'r-', linewidth=2)
527
plt.xlabel('Linear Power')
528
plt.ylabel('Power (dB)')
529
plt.title('Power to dB Conversion')
530
plt.grid(True)
531

532
plt.subplot(2, 3, 3)
533
plt.plot(power_db, power_recovered, 'g-', linewidth=2)
534
plt.xlabel('Power (dB)')
535
plt.ylabel('Recovered Linear Power')
536
plt.title('dB to Power Conversion')
537
plt.grid(True)
538

539
# Magnitude conversions
540
plt.subplot(2, 3, 4)
541
plt.loglog(linear_magnitudes, linear_magnitudes, 'b-', linewidth=2, label='Linear scale')
542
plt.xlabel('Input Magnitude')
543
plt.ylabel('Magnitude')
544
plt.title('Linear Magnitude Scale')
545
plt.grid(True)
546
plt.legend()
547

548
plt.subplot(2, 3, 5)
549
plt.semilogx(linear_magnitudes, magnitude_db, 'r-', linewidth=2)
550
plt.xlabel('Linear Magnitude')
551
plt.ylabel('Magnitude (dB)')
552
plt.title('Magnitude to dB Conversion')
553
plt.grid(True)
554

555
plt.subplot(2, 3, 6)
556
plt.plot(magnitude_db, magnitude_recovered, 'g-', linewidth=2)
557
plt.xlabel('Magnitude (dB)')
558
plt.ylabel('Recovered Linear Magnitude')
559
plt.title('dB to Magnitude Conversion')
560
plt.grid(True)
561

562
plt.tight_layout()
563
plt.show()
564

565
# Demonstrate relationship between power and magnitude dB
566
print("Power vs Magnitude dB Relationship:")
567
print("-" * 40)
568
test_values = [0.1, 1.0, 10.0, 100.0]
569

570
for val in test_values:
571
    pow_db = spectrum.pow2db(val)
572
    mag_db = spectrum.mag2db(np.sqrt(val))  # magnitude = sqrt(power)
573
    print(f"Power: {val:6.1f} -> {pow_db:6.1f} dB")
574
    print(f"  Mag: {np.sqrt(val):6.3f} -> {mag_db:6.1f} dB")
575
    print(f"  Ratio: {pow_db/mag_db:.1f} (should be 2.0)")
576
    print()
577

578
# Practical example: Signal analysis with dB scaling
579
fs = 1000
580
t = np.linspace(0, 1, fs)
581

582
# Multi-tone signal with different amplitudes
583
signal = (2.0 * np.sin(2*np.pi*50*t) +    # Strong component
584
         0.5 * np.sin(2*np.pi*150*t) +    # Medium component  
585
         0.1 * np.sin(2*np.pi*300*t) +    # Weak component
586
         0.02 * np.random.randn(fs))       # Noise
587

588
# Compute PSD and convert to dB
589
f, psd_linear = signal.welch(signal, fs, nperseg=256)
590
psd_db = spectrum.pow2db(psd_linear)
591

592
plt.figure(figsize=(12, 8))
593

594
plt.subplot(2, 2, 1)
595
plt.plot(t[:200], signal[:200])
596
plt.title('Multi-tone Signal')
597
plt.xlabel('Time (s)')
598
plt.ylabel('Amplitude')
599
plt.grid(True)
600

601
plt.subplot(2, 2, 2)
602
plt.plot(f, psd_linear)
603
plt.title('PSD - Linear Scale')
604
plt.xlabel('Frequency (Hz)')
605
plt.ylabel('PSD')
606
plt.grid(True)
607

608
plt.subplot(2, 2, 3)
609
plt.semilogy(f, psd_linear)
610
plt.title('PSD - Log Scale') 
611
plt.xlabel('Frequency (Hz)')
612
plt.ylabel('PSD')
613
plt.grid(True)
614

615
plt.subplot(2, 2, 4)
616
plt.plot(f, psd_db)
617
plt.title('PSD - dB Scale')
618
plt.xlabel('Frequency (Hz)')
619
plt.ylabel('PSD (dB)')
620
plt.grid(True)
621

622
plt.tight_layout()
623
plt.show()
624

625
# Dynamic range comparison
626
linear_range = np.max(psd_linear) / np.min(psd_linear[psd_linear > 0])
627
db_range = np.max(psd_db) - np.min(psd_db[np.isfinite(psd_db)])
628

629
print(f"Dynamic Range Analysis:")
630
print(f"Linear scale: {linear_range:.1e} (ratio)")
631
print(f"dB scale: {db_range:.1f} dB")
632
print(f"Equivalent ratio: {spectrum.db2pow(db_range):.1e}")
633
```
634

635
### Utility Functions and Data Format Conversions
636

637
```python
638
import spectrum
639
import numpy as np
640
import matplotlib.pyplot as plt
641

642
# Demonstrate utility functions
643
print("Utility Functions Demonstration:")
644
print("-" * 35)
645

646
# Next power of 2
647
test_values = [100, 255, 256, 513, 1000, 1024]
648
for val in test_values:
649
    next_pow2 = spectrum.nextpow2(val)
650
    print(f"nextpow2({val:4d}) = {next_pow2:4d} (2^{int(np.log2(next_pow2))})")
651

652
# Generate test PSD for format conversions
653
N = 512
654
freqs = np.linspace(0, 0.5, N//2 + 1)  # One-sided frequencies
655
# Create artificial PSD with multiple peaks
656
psd_onesided = (np.exp(-(freqs-0.1)**2/0.01) + 
657
               0.5*np.exp(-(freqs-0.3)**2/0.005) + 
658
               0.1*np.ones_like(freqs))  # Noise floor
659

660
plt.figure(figsize=(15, 10))
661

662
# Original one-sided PSD
663
plt.subplot(2, 3, 1)
664
plt.plot(freqs, psd_onesided, 'b-', linewidth=2)
665
plt.title('Original One-sided PSD')
666
plt.xlabel('Normalized Frequency')
667
plt.ylabel('PSD')
668
plt.grid(True)
669

670
# Convert to two-sided
671
psd_twosided = spectrum.onesided_2_twosided(psd_onesided)
672
freqs_twosided = np.linspace(-0.5, 0.5, len(psd_twosided))
673

674
plt.subplot(2, 3, 2)
675
plt.plot(freqs_twosided, psd_twosided, 'r-', linewidth=2)
676
plt.title('Two-sided PSD')
677
plt.xlabel('Normalized Frequency')  
678
plt.ylabel('PSD')
679
plt.grid(True)
680

681
# Convert two-sided to center-DC
682
psd_centerdc = spectrum.twosided_2_centerdc(psd_twosided)
683

684
plt.subplot(2, 3, 3)
685
plt.plot(freqs_twosided, psd_centerdc, 'g-', linewidth=2)
686
plt.title('Center-DC Format PSD')
687
plt.xlabel('Normalized Frequency')
688
plt.ylabel('PSD')
689
plt.grid(True)
690

691
# Convert back to two-sided
692
psd_twosided_recovered = spectrum.centerdc_2_twosided(psd_centerdc)
693

694
plt.subplot(2, 3, 4)
695
plt.plot(freqs_twosided, psd_twosided_recovered, 'm--', linewidth=2, label='Recovered')
696
plt.plot(freqs_twosided, psd_twosided, 'r:', alpha=0.7, label='Original')
697
plt.title('Recovered Two-sided PSD')
698
plt.xlabel('Normalized Frequency')
699
plt.ylabel('PSD')
700
plt.legend()
701
plt.grid(True)
702

703
# Convert back to one-sided
704
psd_onesided_recovered = spectrum.twosided_2_onesided(psd_twosided_recovered)
705

706
plt.subplot(2, 3, 5)
707
plt.plot(freqs, psd_onesided_recovered, 'c--', linewidth=2, label='Recovered')
708
plt.plot(freqs, psd_onesided, 'b:', alpha=0.7, label='Original')
709
plt.title('Recovered One-sided PSD')
710
plt.xlabel('Normalized Frequency')
711
plt.ylabel('PSD')
712
plt.legend()
713
plt.grid(True)
714

715
# Circular shift demonstration
716
plt.subplot(2, 3, 6)
717
test_signal = np.sin(2*np.pi*0.1*np.arange(100))
718
shifts = [0, 10, -15, 25]
719
for i, shift in enumerate(shifts):
720
    shifted = spectrum.cshift(test_signal, shift)
721
    plt.plot(shifted[:50], label=f'Shift {shift}', alpha=0.8)
722
plt.title('Circular Shift Examples')
723
plt.xlabel('Sample')
724
plt.ylabel('Amplitude')
725
plt.legend()
726
plt.grid(True)
727

728
plt.tight_layout()
729
plt.show()
730

731
# Verify conversion accuracy
732
print("\nFormat Conversion Accuracy:")
733
print("-" * 30)
734
error_onesided = np.max(np.abs(psd_onesided - psd_onesided_recovered))
735
error_twosided = np.max(np.abs(psd_twosided - psd_twosided_recovered))
736

737
print(f"One-sided round-trip error: {error_onesided:.2e}")
738
print(f"Two-sided round-trip error: {error_twosided:.2e}")
739

740
# Demonstrate FFT shift
741
fft_example = np.fft.fft(np.random.randn(16))
742
print("\nFFT Shift Example:")
743
print("Original FFT order:", np.arange(16))
744
print("After fftshift:   ", np.arange(16)[np.fft.fftshift(np.arange(16)) == np.arange(16)])
745

746
# Show frequency ordering
747
N_fft = 16
748
original_freqs = np.fft.fftfreq(N_fft)
749
shifted_freqs = spectrum.fftshift(original_freqs)
750

751
print(f"Original freq bins: {original_freqs}")
752
print(f"Shifted freq bins:  {shifted_freqs}")
753
```