0
# Linear Prediction
1

2
Linear prediction analysis and parameter conversion utilities for autoregressive modeling and speech processing applications. These functions provide comprehensive conversion between different parameter representations used in digital signal processing and linear predictive coding.
3

4
## Capabilities
5

6
### Levinson-Durbin Algorithm
7

8
Core algorithm for solving the Yule-Walker equations to obtain AR parameters from autocorrelation sequences.
9

10
```python { .api }
11
def LEVINSON(r, order=None, allow_singularity=False):
12
    """
13
    Levinson-Durbin recursion for AR parameter estimation.
14
    
15
    Parameters:
16
    - r: array-like, autocorrelation sequence (r[0] is zero-lag)
17
    - order: int, AR model order (None uses len(r)-1)
18
    - allow_singularity: bool, allow singular autocorrelation matrix
19
    
20
    Returns:
21
    tuple: (AR coefficients array, prediction error float, reflection coefficients array)
22
    """
23

24
def rlevinson(a, efinal):
25
    """
26
    Reverse Levinson-Durbin (AR coefficients to autocorrelation).
27
    
28
    Parameters:
29
    - a: array-like, AR coefficients including a0=1
30
    - efinal: float, final prediction error
31
    
32
    Returns:
33
    array: Autocorrelation sequence
34
    """
35

36
def levdown(anxt, enxt=None):
37
    """
38
    Levinson downward recursion.
39
    
40
    Parameters:
41
    - anxt: array-like, AR coefficients at order n+1
42
    - enxt: float, prediction error at order n+1
43
    
44
    Returns:
45
    tuple: (AR coefficients at order n, prediction error at order n)
46
    """
47

48
def levup(acur, knxt, ecur=None):
49
    """
50
    Levinson upward recursion.
51
    
52
    Parameters:
53
    - acur: array-like, AR coefficients at current order
54
    - knxt: float, next reflection coefficient
55
    - ecur: float, current prediction error
56
    
57
    Returns:
58
    tuple: (AR coefficients at next order, prediction error at next order)
59
    """
60
```
61

62
### Linear Predictive Coding
63

64
High-level LPC analysis function for speech and audio processing applications.
65

66
```python { .api }
67
def lpc(x, N=None):
68
    """
69
    Linear Predictive Coding analysis.
70
    
71
    Parameters:
72
    - x: array-like, input time series data
73
    - N: int, LPC order (None for automatic selection)
74
    
75
    Returns:
76
    tuple: (LPC coefficients array, prediction error float)
77
    """
78
```
79

80
### Autocorrelation Conversions
81

82
Functions for converting between autocorrelation sequences and polynomial coefficients.
83

84
```python { .api }
85
def ac2poly(data):
86
    """
87
    Autocorrelation to polynomial coefficient conversion.
88
    
89
    Parameters:
90
    - data: array-like, autocorrelation sequence
91
    
92
    Returns:
93
    array: Polynomial coefficients
94
    """
95

96
def poly2ac(poly, efinal):
97
    """
98
    Polynomial coefficients to autocorrelation conversion.
99
    
100
    Parameters:
101
    - poly: array-like, polynomial coefficients
102
    - efinal: float, final prediction error
103
    
104
    Returns:
105
    array: Autocorrelation sequence
106
    """
107

108
def ac2rc(data):
109
    """
110
    Autocorrelation to reflection coefficients conversion.
111
    
112
    Parameters:
113
    - data: array-like, autocorrelation sequence
114
    
115
    Returns:
116
    array: Reflection coefficients
117
    """
118

119
def rc2ac(k, R0):
120
    """
121
    Reflection coefficients to autocorrelation conversion.
122
    
123
    Parameters:
124
    - k: array-like, reflection coefficients
125
    - R0: float, zero-lag autocorrelation
126
    
127
    Returns:
128
    array: Autocorrelation sequence
129
    """
130
```
131

132
### Reflection Coefficient Conversions
133

134
Functions for converting between AR parameters and reflection coefficients (PARCOR coefficients).
135

136
```python { .api }
137
def ar2rc(ar):
138
    """
139
    AR parameters to reflection coefficients conversion.
140
    
141
    Parameters:
142
    - ar: array-like, AR coefficients (without leading 1)
143
    
144
    Returns:
145
    array: Reflection coefficients
146
    """
147

148
def rc2poly(kr, r0=None):
149
    """
150
    Reflection coefficients to polynomial coefficients conversion.
151
    
152
    Parameters:
153
    - kr: array-like, reflection coefficients
154
    - r0: float, zero-lag autocorrelation (optional)
155
    
156
    Returns:
157
    array: Polynomial coefficients
158
    """
159

160
def poly2rc(a, efinal):
161
    """
162
    Polynomial coefficients to reflection coefficients conversion.
163
    
164
    Parameters:
165
    - a: array-like, polynomial coefficients
166
    - efinal: float, final prediction error
167
    
168
    Returns:
169
    array: Reflection coefficients
170
    """
171
```
172

173
### Line Spectral Frequencies
174

175
Conversion between polynomial coefficients and line spectral frequencies (LSF), useful for quantization and stability analysis.
176

177
```python { .api }
178
def poly2lsf(a):
179
    """
180
    Polynomial coefficients to line spectral frequencies conversion.
181
    
182
    Parameters:
183
    - a: array-like, polynomial coefficients
184
    
185
    Returns:
186
    array: Line spectral frequencies (normalized to [0, π])
187
    """
188

189
def lsf2poly(lsf):
190
    """
191
    Line spectral frequencies to polynomial coefficients conversion.
192
    
193
    Parameters:
194
    - lsf: array-like, line spectral frequencies
195
    
196
    Returns:
197
    array: Polynomial coefficients
198
    """
199
```
200

201
### Alternative Parameterizations
202

203
Conversions involving inverse sine parameters and log area ratios for robust parameter representation.
204

205
```python { .api }
206
def is2rc(inv_sin):
207
    """
208
    Inverse sine parameters to reflection coefficients conversion.
209
    
210
    Parameters:
211
    - inv_sin: array-like, inverse sine parameters
212
    
213
    Returns:
214
    array: Reflection coefficients
215
    """
216

217
def rc2is(k):
218
    """
219
    Reflection coefficients to inverse sine parameters conversion.
220
    
221
    Parameters:
222
    - k: array-like, reflection coefficients
223
    
224
    Returns:
225
    array: Inverse sine parameters
226
    """
227

228
def rc2lar(k):
229
    """
230
    Reflection coefficients to log area ratios conversion.
231
    
232
    Parameters:
233
    - k: array-like, reflection coefficients
234
    
235
    Returns:
236
    array: Log area ratios
237
    """
238

239
def lar2rc(g):
240
    """
241
    Log area ratios to reflection coefficients conversion.
242
    
243
    Parameters:
244
    - g: array-like, log area ratios
245
    
246
    Returns:
247
    array: Reflection coefficients
248
    """
249
```
250

251
## Usage Examples
252

253
### Basic Linear Prediction Analysis
254

255
```python
256
import spectrum
257
import numpy as np
258
import matplotlib.pyplot as plt
259

260
# Generate AR signal for LPC analysis
261
N = 256
262
# Create AR(3) signal with known coefficients
263
ar_true = np.array([1.0, -2.2137, 2.9403, -1.7720])  # Includes leading 1
264
noise = np.random.randn(N)
265

266
# Generate signal using difference equation
267
signal = np.zeros(N)
268
signal[:3] = noise[:3]  # Initialize
269
for n in range(3, N):
270
    signal[n] = -sum(ar_true[1:] * signal[n-3:n][::-1]) + noise[n]
271

272
# LPC analysis
273
lpc_order = 3
274
lpc_coeffs, pred_error = spectrum.lpc(signal, N=lpc_order)
275

276
print(f"True AR coefficients: {-ar_true[1:]}")
277
print(f"LPC coefficients: {lpc_coeffs[1:]}")
278
print(f"Prediction error: {pred_error}")
279

280
# Compare with autocorrelation method
281
autocorr = spectrum.CORRELATION(signal, maxlags=lpc_order, norm='biased')
282
ar_levinson, err_levinson, rc_levinson = spectrum.LEVINSON(autocorr, order=lpc_order)
283

284
print(f"Levinson AR coeffs: {ar_levinson}")
285
print(f"Levinson error: {err_levinson}")
286
print(f"Reflection coeffs: {rc_levinson}")
287
```
288

289
### Parameter Conversion Chain
290

291
```python
292
import spectrum
293
import numpy as np
294

295
# Start with AR coefficients
296
ar_coeffs = np.array([0.5, -0.3, 0.1])  # AR(3) coefficients
297
print(f"Original AR coefficients: {ar_coeffs}")
298

299
# Convert AR to reflection coefficients
300
rc = spectrum.ar2rc(ar_coeffs)
301
print(f"Reflection coefficients: {rc}")
302

303
# Convert reflection coefficients back to AR
304
ar_reconstructed = spectrum.rc2poly(rc)[1:]  # Remove leading 1
305
print(f"Reconstructed AR: {ar_reconstructed}")
306

307
# Convert to line spectral frequencies
308
poly_with_one = np.concatenate([[1], ar_coeffs])
309
lsf = spectrum.poly2lsf(poly_with_one)
310
print(f"Line spectral frequencies: {lsf}")
311

312
# Convert LSF back to polynomial
313
poly_reconstructed = spectrum.lsf2poly(lsf)
314
ar_from_lsf = poly_reconstructed[1:]
315
print(f"AR from LSF: {ar_from_lsf}")
316

317
# Convert to log area ratios
318
lar = spectrum.rc2lar(rc)
319
print(f"Log area ratios: {lar}")
320

321
# Convert LAR back to reflection coefficients
322
rc_from_lar = spectrum.lar2rc(lar)
323
print(f"RC from LAR: {rc_from_lar}")
324

325
# Verify round-trip conversion accuracy
326
print(f"\nConversion errors:")
327
print(f"AR round-trip error: {np.max(np.abs(ar_coeffs - ar_reconstructed))}")
328
print(f"AR from LSF error: {np.max(np.abs(ar_coeffs - ar_from_lsf))}")
329
print(f"RC from LAR error: {np.max(np.abs(rc - rc_from_lar))}")
330
```
331

332
### Autocorrelation-Based Analysis
333

334
```python
335
import spectrum
336
import numpy as np
337
import matplotlib.pyplot as plt
338

339
# Generate test signal
340
fs = 1000
341
N = 500
342
t = np.arange(N) / fs
343

344
# Multi-component signal
345
signal = (np.sin(2*np.pi*50*t) + 0.5*np.sin(2*np.pi*120*t) + 
346
         0.2*np.random.randn(N))
347

348
# Compute autocorrelation with different methods
349
max_lag = 20
350
autocorr_biased = spectrum.CORRELATION(signal, maxlags=max_lag, norm='biased')
351
autocorr_unbiased = spectrum.CORRELATION(signal, maxlags=max_lag, norm='unbiased')
352

353
# LPC analysis using autocorrelation
354
orders = [5, 10, 15, 20]
355
plt.figure(figsize=(15, 10))
356

357
# Plot signal and autocorrelation
358
plt.subplot(2, 3, 1)
359
plt.plot(t[:200], signal[:200])
360
plt.title('Input Signal')
361
plt.xlabel('Time (s)')
362
plt.ylabel('Amplitude')
363
plt.grid(True)
364

365
plt.subplot(2, 3, 2)
366
lags = np.arange(-max_lag, max_lag+1)
367
plt.plot(lags, autocorr_biased, 'b-', label='Biased')
368
plt.plot(lags, autocorr_unbiased, 'r--', label='Unbiased')
369
plt.title('Autocorrelation Function')
370
plt.xlabel('Lag')
371
plt.ylabel('Correlation')
372
plt.legend()
373
plt.grid(True)
374

375
# LPC analysis for different orders
376
for i, order in enumerate(orders):
377
    plt.subplot(2, 3, 3+i)
378
    
379
    # Use positive lags for Levinson-Durbin
380
    autocorr_pos = autocorr_biased[max_lag:]  # Take positive lags including zero
381
    ar_coeffs, pred_error, rc = spectrum.LEVINSON(autocorr_pos, order=order)
382
    
383
    # Compute and plot frequency response
384
    w = np.linspace(0, np.pi, 512)
385
    # H(w) = 1/A(w) where A(w) is AR polynomial
386
    A_w = np.polyval(np.concatenate([[1], ar_coeffs]), np.exp(1j*w))
387
    H_w = 1 / A_w
388
    psd_ar = np.abs(H_w)**2 * pred_error
389
    
390
    freqs_hz = w * fs / (2*np.pi)
391
    plt.semilogy(freqs_hz, psd_ar, 'b-', linewidth=2)
392
    plt.title(f'LPC Spectrum (Order {order})')
393
    plt.xlabel('Frequency (Hz)')
394
    plt.ylabel('PSD')
395
    plt.xlim(0, fs/2)
396
    plt.grid(True)
397
    
398
    # Print some statistics
399
    print(f"Order {order}: Pred error = {pred_error:.6f}, RC = {rc[:3]}")
400

401
plt.tight_layout()
402
plt.show()
403
```
404

405
### Stability Analysis Using Different Parameterizations
406

407
```python
408
import spectrum
409
import numpy as np
410
import matplotlib.pyplot as plt
411

412
def check_stability(ar_coeffs):
413
    """Check if AR filter is stable (all poles inside unit circle)."""
414
    # Poles are roots of AR polynomial
415
    poly = np.concatenate([[1], ar_coeffs])
416
    poles = np.roots(poly)
417
    return np.all(np.abs(poles) < 1), poles
418

419
# Generate series of AR coefficients with varying stability
420
test_cases = [
421
    np.array([0.5, -0.2]),           # Stable
422
    np.array([0.8, -0.6]),           # Stable  
423
    np.array([1.2, -0.5]),           # Unstable
424
    np.array([-1.8, 0.9]),           # Marginally stable
425
    np.array([0.9, 0.9, -0.5])      # Higher order
426
]
427

428
plt.figure(figsize=(15, 12))
429

430
for i, ar_coeffs in enumerate(test_cases):
431
    # Check stability
432
    is_stable, poles = check_stability(ar_coeffs)
433
    
434
    # Convert to different representations
435
    rc = spectrum.ar2rc(ar_coeffs)
436
    poly_with_one = np.concatenate([[1], ar_coeffs])
437
    
438
    try:
439
        lsf = spectrum.poly2lsf(poly_with_one)
440
        lar = spectrum.rc2lar(rc)
441
        
442
        # Plot pole-zero diagram
443
        plt.subplot(len(test_cases), 4, 4*i + 1)
444
        plt.plot(np.real(poles), np.imag(poles), 'ro', markersize=8)
445
        circle = plt.Circle((0, 0), 1, fill=False, linestyle='--', color='blue')
446
        plt.gca().add_patch(circle)
447
        plt.axis('equal')
448
        plt.xlim(-1.5, 1.5)
449
        plt.ylim(-1.5, 1.5)
450
        plt.title(f'Case {i+1}: {"Stable" if is_stable else "Unstable"}')
451
        plt.xlabel('Real')
452
        plt.ylabel('Imag')
453
        plt.grid(True)
454
        
455
        # Plot reflection coefficients
456
        plt.subplot(len(test_cases), 4, 4*i + 2)
457
        plt.stem(range(len(rc)), rc, basefmt=' ')
458
        plt.axhline(1, color='r', linestyle='--', alpha=0.7)
459
        plt.axhline(-1, color='r', linestyle='--', alpha=0.7)
460
        plt.title('Reflection Coefficients')
461
        plt.ylabel('RC Value')
462
        plt.grid(True)
463
        
464
        # Plot line spectral frequencies
465
        plt.subplot(len(test_cases), 4, 4*i + 3)
466
        plt.stem(range(len(lsf)), lsf/np.pi, basefmt=' ')
467
        plt.title('Line Spectral Frequencies')
468
        plt.ylabel('LSF (normalized)')
469
        plt.ylim(0, 1)
470
        plt.grid(True)
471
        
472
        # Plot log area ratios
473
        plt.subplot(len(test_cases), 4, 4*i + 4)
474
        plt.stem(range(len(lar)), lar, basefmt=' ')
475
        plt.title('Log Area Ratios')
476
        plt.ylabel('LAR')
477
        plt.grid(True)
478
        
479
        # Print parameter ranges for stability analysis
480
        print(f"\nCase {i+1} ({'Stable' if is_stable else 'Unstable'}):")
481
        print(f"  AR coeffs: {ar_coeffs}")
482
        print(f"  RC range: [{np.min(rc):.3f}, {np.max(rc):.3f}]")
483
        print(f"  RC stability: {np.all(np.abs(rc) < 1)}")
484
        print(f"  Poles: {poles}")
485
        
486
    except Exception as e:
487
        print(f"Error processing case {i+1}: {e}")
488

489
plt.tight_layout()
490
plt.show()
491

492
# Demonstrate quantization robustness
493
print("\n" + "="*60)
494
print("Quantization Robustness Comparison")
495
print("="*60)
496

497
# Original AR coefficients  
498
ar_orig = np.array([0.7, -0.4, 0.15])
499
print(f"Original AR: {ar_orig}")
500

501
# Add quantization noise
502
quantization_levels = [8, 6, 4]  # bits
503
for bits in quantization_levels:
504
    # Quantize AR coefficients directly
505
    scale = 2**(bits-1) - 1
506
    ar_quant = np.round(ar_orig * scale) / scale
507
    
508
    # Quantize via reflection coefficients
509
    rc_orig = spectrum.ar2rc(ar_orig)
510
    rc_quant = np.round(rc_orig * scale) / scale
511
    ar_from_rc = spectrum.rc2poly(rc_quant)[1:]
512
    
513
    # Quantize via LSF
514
    poly_orig = np.concatenate([[1], ar_orig])
515
    lsf_orig = spectrum.poly2lsf(poly_orig)
516
    lsf_quant = np.round(lsf_orig * scale) / scale
517
    poly_from_lsf = spectrum.lsf2poly(lsf_quant)
518
    ar_from_lsf = poly_from_lsf[1:]
519
    
520
    # Check stability
521
    stable_direct = check_stability(ar_quant)[0]
522
    stable_rc = check_stability(ar_from_rc)[0]
523
    stable_lsf = check_stability(ar_from_lsf)[0]
524
    
525
    print(f"\n{bits}-bit quantization:")
526
    print(f"  Direct AR:    {ar_quant}, stable: {stable_direct}")
527
    print(f"  Via RC:       {ar_from_rc}, stable: {stable_rc}")
528
    print(f"  Via LSF:      {ar_from_lsf}, stable: {stable_lsf}")
529
```