0
# Mathematical Tools
1

2
Linear algebra utilities, model selection criteria, transfer function analysis, and mathematical tools supporting spectral analysis. These functions provide the mathematical foundation for advanced signal processing and system identification applications.
3

4
## Capabilities
5

6
### Linear Algebra Utilities
7

8
Core linear algebra functions for matrix operations and decompositions used in spectral analysis.
9

10
```python { .api }
11
def CHOLESKY(A, B, method='scipy'):
12
    """
13
    Solve linear system AX=B using Cholesky decomposition.
14
    
15
    Parameters:
16
    - A: array-like, Hermitian positive definite coefficient matrix
17
    - B: array-like, right-hand side matrix or vector
18
    - method: str, computation method ('scipy', 'numpy', 'numpy_solver')
19
    
20
    Returns:
21
    array: Solution matrix X such that AX = B
22
    """
23

24
def pascal(n):
25
    """
26
    Generate Pascal matrix of size n x n.
27
    
28
    Parameters:
29
    - n: int, matrix size
30
    
31
    Returns:
32
    array: Pascal matrix (symmetric positive definite)
33
    """
34

35
def corrmtx(x_input, m, method='autocorrelation'):
36
    """
37
    Generate correlation matrix for linear prediction.
38
    
39
    Parameters:
40
    - x_input: array-like, input data vector
41
    - m: int, matrix order (number of coefficients)
42
    - method: str, matrix construction method ('autocorrelation', 'covariance', 'prewindowed', 'postwindowed')
43
    
44
    Returns:
45
    array: Correlation matrix for linear prediction analysis
46
    """
47

48
def csvd(A):
49
    """
50
    Complex singular value decomposition.
51
    
52
    Parameters:
53
    - A: array-like, input matrix (real or complex)
54
    
55
    Returns:
56
    tuple: (U array, singular values array, V^H array)
57
    """
58
```
59

60
### Model Selection Criteria
61

62
Information-theoretic criteria for automatic model order selection in parametric spectral estimation.
63

64
```python { .api }
65
def AIC(N, rho, k):
66
    """
67
    Akaike Information Criterion for model order selection.
68
    
69
    Parameters:
70
    - N: int, data length
71
    - rho: float, prediction error (noise variance)
72
    - k: int, model order (number of parameters)
73
    
74
    Returns:
75
    float: AIC value (lower is better)
76
    """
77

78
def AICc(N, rho, k, norm=True):
79
    """
80
    Corrected AIC for small sample sizes.
81
    
82
    Parameters:
83
    - N: int, data length
84
    - rho: float, prediction error
85
    - k: int, model order
86
    - norm: bool, apply normalization factor
87
    
88
    Returns:
89
    float: Corrected AIC value
90
    """
91

92
def MDL(N, rho, k):
93
    """
94
    Minimum Description Length criterion.
95
    
96
    Parameters:
97
    - N: int, data length
98
    - rho: float, prediction error
99
    - k: int, model order
100
    
101
    Returns:
102
    float: MDL value (lower is better)
103
    """
104

105
def FPE(N, rho, k=None):
106
    """
107
    Final Prediction Error criterion.
108
    
109
    Parameters:
110
    - N: int, data length
111
    - rho: float, prediction error
112
    - k: int, model order (None for simple FPE)
113
    
114
    Returns:
115
    float: FPE value
116
    """
117

118
def KIC(N, rho, k):
119
    """
120
    Kullback Information Criterion.
121
    
122
    Parameters:
123
    - N: int, data length
124
    - rho: float, prediction error
125
    - k: int, model order
126
    
127
    Returns:
128
    float: KIC value
129
    """
130

131
def AKICc(N, rho, k):
132
    """
133
    Asymmetric KIC corrected for small samples.
134
    
135
    Parameters:
136
    - N: int, data length
137
    - rho: float, prediction error
138
    - k: int, model order
139
    
140
    Returns:
141
    float: AKICc value
142
    """
143

144
def CAT(N, rho, k):
145
    """
146
    Criterion Autoregressive Transfer function.
147
    
148
    Parameters:
149
    - N: int, data length
150
    - rho: float, prediction error
151
    - k: int, model order
152
    
153
    Returns:
154
    float: CAT value
155
    """
156
```
157

158
### Automatic Model Order Selection
159

160
Class-based interface for automatic model order selection using various criteria.
161

162
```python { .api }
163
class Criteria(name, N):
164
    """
165
    Automatic order selection for parametric methods.
166
    
167
    Parameters:
168
    - name: str, criterion name ('AIC', 'MDL', 'FPE', 'CAT', etc.)
169
    - N: int, data length
170
    
171
    Methods:
172
    - __call__(rho, k): compute criterion value
173
    - find_order(rho_values, orders): find optimal order
174
    """
175
```
176

177
### Transfer Function Analysis
178

179
Functions for converting between different transfer function representations and performing system analysis.
180

181
```python { .api }
182
def tf2zp(b, a):
183
    """
184
    Transfer function to zeros and poles conversion.
185
    
186
    Parameters:
187
    - b: array-like, numerator polynomial coefficients
188
    - a: array-like, denominator polynomial coefficients
189
    
190
    Returns:
191
    tuple: (zeros array, poles array)
192
    """
193

194
def zp2tf(z, p, k):
195
    """
196
    Zeros and poles to transfer function conversion.
197
    
198
    Parameters:
199
    - z: array-like, zeros of the transfer function
200
    - p: array-like, poles of the transfer function  
201
    - k: float, system gain
202
    
203
    Returns:
204
    tuple: (numerator coefficients array, denominator coefficients array)
205
    """
206

207
def tf2zpk(b, a):
208
    """
209
    Transfer function to zeros, poles, and gain conversion.
210
    
211
    Parameters:
212
    - b: array-like, numerator polynomial coefficients
213
    - a: array-like, denominator polynomial coefficients
214
    
215
    Returns:
216
    tuple: (zeros array, poles array, gain float)
217
    """
218

219
def zpk2tf(z, p, k):
220
    """
221
    Zeros, poles, and gain to transfer function conversion.
222
    
223
    Parameters:
224
    - z: array-like, zeros of the transfer function
225
    - p: array-like, poles of the transfer function
226
    - k: float, system gain
227
    
228
    Returns:
229
    tuple: (numerator coefficients array, denominator coefficients array)
230
    """
231

232
def eqtflength(b, a):
233
    """
234
    Equalize lengths of transfer function numerator and denominator.
235
    
236
    Parameters:
237
    - b: array-like, numerator coefficients
238
    - a: array-like, denominator coefficients
239
    
240
    Returns:
241
    tuple: (equalized numerator array, equalized denominator array)
242
    """
243
```
244

245
### Advanced Transfer Function Conversions
246

247
Functions for state-space and lattice filter representations.
248

249
```python { .api }
250
def zpk2ss(z, p, k):
251
    """
252
    Convert zeros, poles, gain to state-space representation.
253
    
254
    Parameters:
255
    - z: array-like, zeros of the transfer function
256
    - p: array-like, poles of the transfer function
257
    - k: float, system gain
258
    
259
    Returns:
260
    tuple: (A matrix, B vector, C vector, D scalar) state-space matrices
261
    """
262

263
def ss2zpk(a, b, c, d, input=0):
264
    """
265
    Convert state-space to zeros, poles, gain representation.
266
    
267
    Parameters:
268
    - a: array-like, state matrix
269
    - b: array-like, input matrix
270
    - c: array-like, output matrix
271
    - d: array-like, feedthrough matrix
272
    - input: int, input channel for MIMO systems
273
    
274
    Returns:
275
    tuple: (zeros array, poles array, gain float)
276
    """
277

278
def tf2latc(num=[1.], den=[1.]):
279
    """
280
    Convert transfer function to lattice filter coefficients.
281
    
282
    Parameters:
283
    - num: array-like, numerator coefficients (default [1.])
284
    - den: array-like, denominator coefficients (default [1.])
285
    
286
    Returns:
287
    array: Lattice reflection coefficients
288
    """
289

290
def allpole2latc(num, den):
291
    """
292
    Convert all-pole transfer function to lattice form.
293
    
294
    Parameters:
295
    - num: array-like, numerator coefficients
296
    - den: array-like, denominator coefficients
297
    
298
    Returns:
299
    array: Lattice coefficients for all-pole system
300
    """
301
```
302

303
### Toeplitz System Solver
304

305
Specialized solver for Hermitian Toeplitz systems commonly encountered in signal processing.
306

307
```python { .api }
308
def HERMTOEP(T0, T, Z):
309
    """
310
    Solve Hermitian Toeplitz system efficiently.
311
    
312
    Parameters:
313
    - T0: complex, Toeplitz matrix diagonal element
314
    - T: array-like, first row/column of Toeplitz matrix (excluding T0)
315
    - Z: array-like, right-hand side vector
316
    
317
    Returns:
318
    array: Solution vector for Toeplitz system
319
    """
320
```
321

322
## Usage Examples
323

324
### Model Order Selection
325

326
```python
327
import spectrum
328
import numpy as np
329
import matplotlib.pyplot as plt
330

331
# Generate AR signal with known order
332
np.random.seed(42)
333
N = 200
334
true_order = 4
335

336
# AR coefficients for a 4th-order system
337
ar_true = np.array([1.0, -2.7607, 3.8106, -2.6535, 0.9238])
338

339
# Generate AR signal
340
noise = np.random.randn(N)
341
signal = np.zeros(N)
342
signal[:true_order] = noise[:true_order]
343

344
for n in range(true_order, N):
345
    signal[n] = -np.sum(ar_true[1:] * signal[n-true_order:n][::-1]) + noise[n]
346

347
# Test different model orders
348
max_order = 15
349
orders = range(1, max_order + 1)
350
criteria_values = {
351
    'AIC': [],
352
    'MDL': [],
353
    'FPE': [],
354
    'CAT': [],
355
    'AICc': []
356
}
357

358
# Compute AR parameters and criteria for each order
359
for order in orders:
360
    # Estimate AR parameters using Yule-Walker
361
    ar_coeffs, rho, _ = spectrum.aryule(signal, order)
362
    
363
    # Compute various criteria
364
    criteria_values['AIC'].append(spectrum.AIC(N, rho, order))
365
    criteria_values['MDL'].append(spectrum.MDL(N, rho, order))
366
    criteria_values['FPE'].append(spectrum.FPE(N, rho, order))
367
    criteria_values['CAT'].append(spectrum.CAT(N, rho, order))
368
    criteria_values['AICc'].append(spectrum.AICc(N, rho, order))
369

370
# Find optimal orders
371
optimal_orders = {}
372
for name, values in criteria_values.items():
373
    optimal_orders[name] = orders[np.argmin(values)]
374

375
# Plot results
376
plt.figure(figsize=(15, 10))
377

378
plt.subplot(2, 3, 1)
379
plt.plot(signal)
380
plt.title(f'Generated AR({true_order}) Signal')
381
plt.xlabel('Sample')
382
plt.ylabel('Amplitude')
383
plt.grid(True)
384

385
# Plot criteria
386
for i, (name, values) in enumerate(criteria_values.items()):
387
    plt.subplot(2, 3, 2 + i)
388
    plt.plot(orders, values, 'o-', linewidth=2, markersize=6)
389
    plt.axvline(true_order, color='r', linestyle='--', alpha=0.7, label=f'True order = {true_order}')
390
    plt.axvline(optimal_orders[name], color='g', linestyle=':', alpha=0.7, 
391
               label=f'Optimal = {optimal_orders[name]}')
392
    plt.title(f'{name} Criterion')
393
    plt.xlabel('Model Order')
394
    plt.ylabel(f'{name} Value')
395
    plt.legend()
396
    plt.grid(True)
397

398
plt.tight_layout()
399
plt.show()
400

401
# Print results summary
402
print("Model Order Selection Results:")
403
print("-" * 40)
404
print(f"True AR order: {true_order}")
405
print("\nCriterion    Optimal Order")
406
print("-" * 25)
407
for name, optimal in optimal_orders.items():
408
    status = "✓" if optimal == true_order else "✗"
409
    print(f"{name:8s}     {optimal:2d}          {status}")
410

411
# Demonstrate Criteria class usage
412
print("\nUsing Criteria class:")
413
print("-" * 25)
414
aic_criterion = spectrum.Criteria('AIC', N)
415
for order in [3, 4, 5]:
416
    ar_coeffs, rho, _ = spectrum.aryule(signal, order)
417
    aic_val = aic_criterion(rho, order)  # This should work if implemented
418
    print(f"Order {order}: AIC = {aic_val:.3f}")
419
```
420

421
### Linear Algebra and Matrix Operations
422

423
```python
424
import spectrum
425
import numpy as np
426
import matplotlib.pyplot as plt
427

428
# Demonstrate Cholesky solver for linear systems
429
print("Cholesky Decomposition Example:")
430
print("-" * 35)
431

432
# Create a positive definite matrix (correlation matrix)
433
N = 100
434
n = np.arange(N)
435
signal = np.sin(2*np.pi*0.1*n) + 0.1*np.random.randn(N)
436

437
# Create autocorrelation matrix using different methods
438
methods = ['autocorrelation', 'covariance']
439
order = 10
440

441
plt.figure(figsize=(15, 8))
442

443
for i, method in enumerate(methods):
444
    # Generate correlation matrix
445
    R = spectrum.corrmtx(signal, order-1, method=method)
446
    
447
    plt.subplot(2, 4, 1 + i)
448
    plt.imshow(np.real(R), cmap='coolwarm')
449
    plt.colorbar()
450
    plt.title(f'Correlation Matrix ({method})')
451
    
452
    plt.subplot(2, 4, 3 + i)
453
    plt.imshow(np.imag(R), cmap='coolwarm')
454
    plt.colorbar()
455
    plt.title(f'Imaginary Part ({method})')
456
    
457
    # Test Cholesky solver
458
    b = np.random.randn(R.shape[0])  # Random right-hand side
459
    
460
    try:
461
        # Solve using Cholesky decomposition
462
        x_chol = spectrum.CHOLESKY(R, b, method='scipy')
463
        
464
        # Solve using standard method for comparison
465
        x_std = np.linalg.solve(R, b)
466
        
467
        error = np.linalg.norm(x_chol - x_std)
468
        print(f"{method} method:")
469
        print(f"  Matrix size: {R.shape}")
470
        print(f"  Condition number: {np.linalg.cond(R):.2e}")
471
        print(f"  Cholesky error: {error:.2e}")
472
        print(f"  Matrix is positive definite: {np.all(np.linalg.eigvals(R) > 0)}")
473
        
474
    except Exception as e:
475
        print(f"Cholesky failed for {method}: {e}")
476

477
# Pascal matrix example
478
plt.subplot(2, 4, 5)
479
pascal_mat = spectrum.pascal(8)
480
plt.imshow(pascal_mat, cmap='Blues')
481
plt.colorbar()
482
plt.title('Pascal Matrix')
483

484
# Complex SVD demonstration
485
plt.subplot(2, 4, 6)
486
# Create complex test matrix
487
A_complex = np.random.randn(6, 4) + 1j * np.random.randn(6, 4)
488
U, s, Vh = spectrum.csvd(A_complex)
489

490
plt.semilogy(s, 'o-', linewidth=2, markersize=8)
491
plt.title('Singular Values (Complex SVD)')
492
plt.xlabel('Index')
493
plt.ylabel('Singular Value')
494
plt.grid(True)
495

496
plt.tight_layout()
497
plt.show()
498

499
# Verify SVD reconstruction
500
A_reconstructed = U @ np.diag(s) @ Vh
501
reconstruction_error = np.linalg.norm(A_complex - A_reconstructed)
502
print(f"\nComplex SVD verification:")
503
print(f"Original matrix shape: {A_complex.shape}")
504
print(f"Reconstruction error: {reconstruction_error:.2e}")
505
```
506

507
### Transfer Function Analysis
508

509
```python
510
import spectrum
511
import numpy as np
512
import matplotlib.pyplot as plt
513

514
# Define a test transfer function (2nd order low-pass filter)
515
# H(s) = 1 / (s^2 + 2*ζ*ωn*s + ωn^2)
516
wn = 2*np.pi*10  # Natural frequency = 10 Hz
517
zeta = 0.7       # Damping ratio
518

519
# Continuous-time coefficients
520
num_s = [wn**2]
521
den_s = [1, 2*zeta*wn, wn**2]
522

523
# Convert to discrete-time using bilinear transform (approximation)
524
fs = 100  # Sampling frequency
525
T = 1/fs
526

527
# Bilinear transform: s = 2/T * (z-1)/(z+1)
528
# This is a simplified example - normally use scipy.signal.bilinear
529
num_z = num_s
530
den_z = den_s  # Simplified for demonstration
531

532
# Convert to zeros, poles, gain representation
533
zeros, poles, gain = spectrum.tf2zpk(num_z, den_z)
534

535
print("Transfer Function Analysis:")
536
print("-" * 30)
537
print(f"Numerator coefficients: {num_z}")
538
print(f"Denominator coefficients: {den_z}")
539
print(f"Zeros: {zeros}")
540
print(f"Poles: {poles}")
541
print(f"Gain: {gain}")
542

543
# Convert back to transfer function
544
num_recovered, den_recovered = spectrum.zpk2tf(zeros, poles, gain)
545
print(f"Recovered numerator: {num_recovered}")
546
print(f"Recovered denominator: {den_recovered}")
547

548
# Verify conversion accuracy
549
num_error = np.max(np.abs(np.array(num_z) - num_recovered[:len(num_z)]))
550
den_error = np.max(np.abs(np.array(den_z) - den_recovered[:len(den_z)]))
551
print(f"Numerator error: {num_error:.2e}")
552
print(f"Denominator error: {den_error:.2e}")
553

554
# Pole-zero plot and frequency response
555
plt.figure(figsize=(15, 8))
556

557
# Pole-zero plot
558
plt.subplot(2, 3, 1)
559
plt.plot(np.real(zeros), np.imag(zeros), 'bo', markersize=10, label='Zeros')
560
plt.plot(np.real(poles), np.imag(poles), 'rx', markersize=10, label='Poles')
561

562
# Unit circle for discrete-time systems
563
theta = np.linspace(0, 2*np.pi, 100)
564
plt.plot(np.cos(theta), np.sin(theta), 'k--', alpha=0.5, label='Unit Circle')
565

566
plt.axis('equal')
567
plt.xlabel('Real Part')
568
plt.ylabel('Imaginary Part')
569
plt.title('Pole-Zero Plot')
570
plt.legend()
571
plt.grid(True)
572

573
# Frequency response
574
w = np.linspace(0, np.pi, 512)  # Normalized frequency
575
H = np.polyval(num_recovered, np.exp(1j*w)) / np.polyval(den_recovered, np.exp(1j*w))
576

577
plt.subplot(2, 3, 2)
578
plt.plot(w/(2*np.pi)*fs, 20*np.log10(np.abs(H)))
579
plt.xlabel('Frequency (Hz)')
580
plt.ylabel('Magnitude (dB)')
581
plt.title('Frequency Response - Magnitude')
582
plt.grid(True)
583

584
plt.subplot(2, 3, 3)
585
plt.plot(w/(2*np.pi)*fs, np.angle(H)*180/np.pi)
586
plt.xlabel('Frequency (Hz)')
587
plt.ylabel('Phase (degrees)')
588
plt.title('Frequency Response - Phase')
589
plt.grid(True)
590

591
# Test lattice conversion for all-pole system (AR filter)
592
ar_coeffs = [1, -1.5, 0.7]  # AR polynomial
593
lattice_coeffs = spectrum.tf2latc(num=[1], den=ar_coeffs)
594

595
plt.subplot(2, 3, 4)
596
plt.stem(range(len(lattice_coeffs)), lattice_coeffs, basefmt=' ')
597
plt.title('Lattice Coefficients')
598
plt.xlabel('Stage')
599
plt.ylabel('Reflection Coefficient')
600
plt.grid(True)
601

602
# System stability analysis
603
plt.subplot(2, 3, 5)
604
# Plot poles in z-plane for stability analysis
605
pole_magnitudes = np.abs(poles)
606
stable = np.all(pole_magnitudes < 1)
607

608
colors = ['green' if mag < 1 else 'red' for mag in pole_magnitudes]
609
plt.scatter(np.real(poles), np.imag(poles), c=colors, s=100, alpha=0.7)
610
plt.plot(np.cos(theta), np.sin(theta), 'k--', alpha=0.5)
611
plt.axis('equal')
612
plt.xlim(-1.5, 1.5)
613
plt.ylim(-1.5, 1.5)
614
plt.xlabel('Real Part')
615
plt.ylabel('Imaginary Part')
616
plt.title(f'System Stability ({"Stable" if stable else "Unstable"})')
617
plt.grid(True)
618

619
# Group delay
620
plt.subplot(2, 3, 6)
621
# Approximate group delay calculation
622
w_gd = w[1:-1]
623
phase = np.angle(H)
624
group_delay = -np.diff(phase) / np.diff(w)
625
group_delay = np.append(group_delay, group_delay[-1])  # Extend for plotting
626

627
plt.plot(w_gd/(2*np.pi)*fs, group_delay[:-1])
628
plt.xlabel('Frequency (Hz)')
629
plt.ylabel('Group Delay (samples)')
630
plt.title('Group Delay')
631
plt.grid(True)
632

633
plt.tight_layout()
634
plt.show()
635

636
# Print stability analysis
637
print(f"\nStability Analysis:")
638
print(f"All poles inside unit circle: {stable}")
639
for i, (pole, mag) in enumerate(zip(poles, pole_magnitudes)):
640
    status = "Stable" if mag < 1 else "Unstable"
641
    print(f"Pole {i+1}: {pole:.3f}, |pole| = {mag:.3f} ({status})")
642
```
643

644
### Advanced Mathematical Tools
645

646
```python
647
import spectrum
648
import numpy as np
649
import matplotlib.pyplot as plt
650

651
# Demonstrate Toeplitz solver
652
print("Toeplitz System Solver:")
653
print("-" * 25)
654

655
# Create Hermitian Toeplitz matrix
656
n = 8
657
T0 = 2.0  # Diagonal element
658
T = np.array([1.5, 0.8, 0.3, 0.1, 0.05, 0.02, 0.01])  # First row/column
659

660
# Construct full Toeplitz matrix for verification
661
from scipy.linalg import toeplitz
662
full_T = toeplitz([T0] + list(T), [T0] + list(T))
663

664
# Random right-hand side
665
Z = np.random.randn(n) + 1j * np.random.randn(n)
666

667
# Solve using specialized Toeplitz solver
668
try:
669
    X_toep = spectrum.HERMTOEP(T0, T, Z)
670
    print(f"Toeplitz solver successful")
671
    
672
    # Verify solution
673
    residual = np.linalg.norm(full_T @ X_toep - Z)
674
    print(f"Residual error: {residual:.2e}")
675
    
676
except Exception as e:
677
    print(f"Toeplitz solver failed: {e}")
678
    # Fall back to standard solver
679
    X_toep = np.linalg.solve(full_T, Z)
680

681
# Compare computational complexity demonstration
682
sizes = [8, 16, 32, 64, 128]
683
times_standard = []
684
times_toeplitz = []
685

686
import time
687

688
plt.figure(figsize=(12, 8))
689

690
plt.subplot(2, 2, 1)
691
plt.imshow(np.real(full_T), cmap='coolwarm')
692
plt.colorbar()
693
plt.title('Toeplitz Matrix (Real Part)')
694

695
plt.subplot(2, 2, 2)
696
solution_error = np.abs(X_toep - np.linalg.solve(full_T, Z))
697
plt.semilogy(solution_error, 'o-')
698
plt.title('Solution Error (Toeplitz vs Standard)')
699
plt.xlabel('Element Index')
700
plt.ylabel('Absolute Error')
701
plt.grid(True)
702

703
# Model selection criteria comparison on real data
704
# Generate test signals with different characteristics
705
plt.subplot(2, 2, 3)
706

707
# Generate AR signals with different orders
708
true_orders = [2, 4, 6, 8]
709
N = 150
710
test_orders = range(1, 12)
711

712
criteria_comparison = {}
713
for criterion in ['AIC', 'MDL', 'FPE']:
714
    criteria_comparison[criterion] = []
715

716
for true_order in true_orders:
717
    # Generate random AR coefficients
718
    ar_coeffs = np.random.randn(true_order + 1)
719
    ar_coeffs[0] = 1  # Normalize
720
    
721
    # Ensure stability by scaling if needed
722
    poles = np.roots(ar_coeffs)
723
    if np.any(np.abs(poles) >= 1):
724
        ar_coeffs[1:] *= 0.9 / np.max(np.abs(poles))
725
    
726
    # Generate signal
727
    noise = np.random.randn(N)
728
    signal = np.zeros(N)
729
    signal[:true_order] = noise[:true_order]
730
    
731
    for n in range(true_order, N):
732
        signal[n] = -np.sum(ar_coeffs[1:] * signal[n-true_order:n][::-1]) + noise[n]
733
    
734
    # Find optimal orders using different criteria
735
    for criterion in criteria_comparison:
736
        criterion_values = []
737
        for order in test_orders:
738
            _, rho, _ = spectrum.aryule(signal, order)
739
            if criterion == 'AIC':
740
                val = spectrum.AIC(N, rho, order)
741
            elif criterion == 'MDL':
742
                val = spectrum.MDL(N, rho, order)
743
            elif criterion == 'FPE':
744
                val = spectrum.FPE(N, rho, order)
745
            criterion_values.append(val)
746
        
747
        optimal_order = test_orders[np.argmin(criterion_values)]
748
        criteria_comparison[criterion].append(optimal_order)
749

750
# Plot criterion performance
751
x_pos = np.arange(len(true_orders))
752
width = 0.25
753

754
for i, criterion in enumerate(criteria_comparison):
755
    plt.bar(x_pos + i*width, criteria_comparison[criterion], width, 
756
           label=criterion, alpha=0.7)
757

758
plt.bar(x_pos + len(criteria_comparison)*width, true_orders, width, 
759
       label='True Order', alpha=0.7, color='red')
760

761
plt.xlabel('Test Signal')
762
plt.ylabel('Selected Order')
763
plt.title('Criterion Performance Comparison')
764
plt.xticks(x_pos + width, [f'AR({order})' for order in true_orders])
765
plt.legend()
766
plt.grid(True, alpha=0.3)
767

768
# Statistical analysis of criterion performance
769
plt.subplot(2, 2, 4)
770
accuracy_scores = {}
771
for criterion in criteria_comparison:
772
    correct_selections = sum(1 for pred, true in zip(criteria_comparison[criterion], true_orders) 
773
                           if pred == true)
774
    accuracy_scores[criterion] = correct_selections / len(true_orders) * 100
775

776
criteria_names = list(accuracy_scores.keys())
777
accuracy_values = list(accuracy_scores.values())
778

779
plt.bar(criteria_names, accuracy_values, color=['blue', 'green', 'orange'])
780
plt.ylabel('Accuracy (%)')
781
plt.title('Order Selection Accuracy')
782
plt.ylim(0, 100)
783
for i, v in enumerate(accuracy_values):
784
    plt.text(i, v + 2, f'{v:.1f}%', ha='center')
785
plt.grid(True, alpha=0.3)
786

787
plt.tight_layout()
788
plt.show()
789

790
print(f"\nCriterion Performance Summary:")
791
print(f"-" * 35)
792
for criterion, accuracy in accuracy_scores.items():
793
    print(f"{criterion:8s}: {accuracy:5.1f}% correct selections")
794

795
# Demonstrate eigenvalue computations for correlation matrices
796
print(f"\nCorrelation Matrix Properties:")
797
print(f"-" * 32)
798

799
# Generate signal and correlation matrix
800
test_signal = np.random.randn(100)
801
R = spectrum.corrmtx(test_signal, 10, method='autocorrelation')
802

803
eigenvals = np.linalg.eigvals(R)
804
condition_number = np.max(eigenvals) / np.min(eigenvals)
805
trace = np.trace(R)
806
determinant = np.linalg.det(R)
807

808
print(f"Matrix size: {R.shape}")
809
print(f"Condition number: {condition_number:.2e}")
810
print(f"Trace: {trace:.3f}")
811
print(f"Determinant: {determinant:.3e}")
812
print(f"Positive definite: {np.all(eigenvals > 0)}")
813
print(f"Eigenvalue range: [{np.min(eigenvals):.3e}, {np.max(eigenvals):.3e}]")
814
```