0
# Eigenvalue Methods
1

2
High-resolution spectral estimation using eigenvalue decomposition of the data covariance matrix. These methods excel at resolving closely spaced frequency components and provide superior performance for short data records, particularly in the presence of noise.
3

4
## Capabilities
5

6
### MUSIC Method
7

8
Multiple Signal Classification (MUSIC) method provides high-resolution pseudospectrum estimation by exploiting the orthogonality between signal and noise subspaces.
9

10
```python { .api }
11
def music(X, IP, NSIG=None, NFFT=4096, threshold=None, criteria='aic', verbose=False):
12
    """
13
    MUSIC pseudospectrum estimation.
14
    
15
    Parameters:
16
    - X: array-like, input time series data
17
    - IP: int, total number of sinusoids + noise sources
18
    - NSIG: int, number of sinusoidal sources (auto-detected if None)
19
    - NFFT: int, FFT length for pseudospectrum computation
20
    - threshold: float, threshold for automatic signal detection
21
    - criteria: str, model order selection criteria ('aic', 'mdl')
22
    - verbose: bool, enable verbose output
23
    
24
    Returns:
25
    array: MUSIC pseudospectrum values
26
    """
27

28
class pmusic(data, IP, NSIG=None, NFFT=None, sampling=1., scale_by_freq=False):
29
    """
30
    MUSIC pseudospectrum class.
31
    
32
    Parameters:
33
    - data: array-like, input time series
34
    - IP: int, total number of complex sinusoids + noise
35
    - NSIG: int, number of complex sinusoids (auto if None)
36
    - NFFT: int, FFT length for pseudospectrum computation
37
    - sampling: float, sampling frequency
38
    - scale_by_freq: bool, scale pseudospectrum by frequency
39
    """
40
```
41

42
### Eigenvector Method
43

44
Eigenvector method provides an alternative eigenvalue-based approach with different weighting of the noise subspace eigenvectors.
45

46
```python { .api }
47
def ev(X, IP, NSIG=None, NFFT=4096, threshold=None, criteria='aic', verbose=False):
48
    """
49
    Eigenvector method pseudospectrum estimation.
50
    
51
    Parameters:
52
    - X: array-like, input time series data
53
    - IP: int, total number of sinusoids + noise sources
54
    - NSIG: int, number of sinusoidal sources (auto-detected if None)
55
    - NFFT: int, FFT length for pseudospectrum computation
56
    - threshold: float, threshold for automatic signal detection
57
    - criteria: str, model order selection criteria ('aic', 'mdl')
58
    - verbose: bool, enable verbose output
59
    
60
    Returns:
61
    array: Eigenvector pseudospectrum values
62
    """
63

64
class pev(data, IP, NSIG=None, NFFT=None, sampling=1., scale_by_freq=False):
65
    """
66
    Eigenvector pseudospectrum class.
67
    
68
    Parameters:
69
    - data: array-like, input time series
70
    - IP: int, total number of complex sinusoids + noise  
71
    - NSIG: int, number of complex sinusoids (auto if None)
72
    - NFFT: int, FFT length for pseudospectrum computation
73
    - sampling: float, sampling frequency
74
    - scale_by_freq: bool, scale pseudospectrum by frequency
75
    """
76
```
77

78
### Minimum Variance Method
79

80
Minimum variance distortionless response method provides adaptive beamforming-based spectral estimation.
81

82
```python { .api }
83
def minvar(X, order, sampling=1., NFFT=4096):
84
    """
85
    Minimum variance pseudospectrum estimation.
86
    
87
    Parameters:
88
    - X: array-like, input time series data
89
    - order: int, order of the minimum variance filter
90
    - sampling: float, sampling frequency
91
    - NFFT: int, FFT length for pseudospectrum computation
92
    
93
    Returns:
94
    array: Minimum variance pseudospectrum values
95
    """
96

97
class pminvar(data, order, NFFT=None, sampling=1., scale_by_freq=False):
98
    """
99
    Minimum variance pseudospectrum class.
100
    
101
    Parameters:
102
    - data: array-like, input time series
103
    - order: int, order of the minimum variance filter
104
    - NFFT: int, FFT length for pseudospectrum computation
105
    - sampling: float, sampling frequency
106
    - scale_by_freq: bool, scale pseudospectrum by frequency
107
    """
108
```
109

110
### General Eigenvalue Method
111

112
Unified interface for eigenvalue-based methods allowing selection between MUSIC and eigenvector methods.
113

114
```python { .api }
115
def eigen(X, P, NSIG=None, method='music', threshold=None, NFFT=4096, criteria='aic', verbose=False):
116
    """
117
    General eigenvalue-based pseudospectrum estimation.
118
    
119
    Parameters:
120
    - X: array-like, input time series data
121
    - P: int, total number of sinusoids + noise sources
122
    - NSIG: int, number of sinusoidal sources (auto if None)
123
    - method: str, method selection ('music', 'ev')
124
    - threshold: float, threshold for signal detection
125
    - NFFT: int, FFT length for pseudospectrum
126
    - criteria: str, order selection criteria ('aic', 'mdl')
127
    - verbose: bool, enable verbose output
128
    
129
    Returns:
130
    array: Pseudospectrum values
131
    """
132
```
133

134
### Model Order Selection
135

136
Automatic determination of signal and noise subspace dimensions using information-theoretic criteria.
137

138
```python { .api }
139
def aic_eigen(s, N):
140
    """
141
    AIC criterion for eigenvalue methods.
142
    
143
    Parameters:
144
    - s: array-like, eigenvalues in descending order
145
    - N: int, data length
146
    
147
    Returns:
148
    array: AIC values for different model orders
149
    """
150

151
def mdl_eigen(s, N):
152
    """
153
    MDL criterion for eigenvalue methods.
154
    
155
    Parameters:
156
    - s: array-like, eigenvalues in descending order
157
    - N: int, data length
158
    
159
    Returns:
160
    array: MDL values for different model orders
161
    """
162
```
163

164
## Usage Examples
165

166
### MUSIC Spectral Analysis
167

168
```python
169
import spectrum
170
import numpy as np
171
import matplotlib.pyplot as plt
172

173
# Generate signal with two closely spaced sinusoids in noise
174
N = 64
175
n = np.arange(N)
176
f1, f2 = 0.25, 0.27  # Closely spaced normalized frequencies
177
SNR_dB = 10
178

179
# Create signal
180
signal = (np.exp(1j * 2 * np.pi * f1 * n) + 
181
         np.exp(1j * 2 * np.pi * f2 * n))
182
noise_power = 10**(-SNR_dB/10) * np.mean(np.abs(signal)**2)
183
noise = np.sqrt(noise_power/2) * (np.random.randn(N) + 1j*np.random.randn(N))
184
noisy_signal = signal + noise
185

186
# MUSIC analysis
187
IP = 15  # Model order (number of complex exponentials + noise)
188
NSIG = 2  # Number of signal components
189

190
# Automatic signal number detection
191
p_music_auto = spectrum.pmusic(noisy_signal, IP=IP)
192
print(f"Auto-detected signals: {p_music_auto.NSIG}")
193

194
# Manual signal specification  
195
p_music = spectrum.pmusic(noisy_signal, IP=IP, NSIG=NSIG)
196

197
# Compare with eigenvector method
198
p_ev = spectrum.pev(noisy_signal, IP=IP, NSIG=NSIG)
199

200
# Plot results
201
freqs = np.linspace(0, 1, len(p_music.psd))
202
plt.figure(figsize=(12, 8))
203

204
plt.subplot(2, 1, 1)
205
plt.plot(n, np.real(noisy_signal), 'b-', alpha=0.7, label='Real part')
206
plt.plot(n, np.imag(noisy_signal), 'r-', alpha=0.7, label='Imag part')
207
plt.title(f'Noisy Signal: Two Complex Sinusoids (f1={f1}, f2={f2})')
208
plt.xlabel('Sample')
209
plt.ylabel('Amplitude')
210
plt.legend()
211
plt.grid(True)
212

213
plt.subplot(2, 1, 2)
214
plt.plot(freqs, 10*np.log10(p_music.psd), 'b-', linewidth=2, label='MUSIC')
215
plt.plot(freqs, 10*np.log10(p_ev.psd), 'r--', linewidth=2, label='Eigenvector')
216
plt.axvline(f1, color='k', linestyle=':', alpha=0.7, label='True f1')
217
plt.axvline(f2, color='k', linestyle=':', alpha=0.7, label='True f2')
218
plt.xlabel('Normalized Frequency')
219
plt.ylabel('Pseudospectrum (dB)')
220
plt.title('MUSIC vs Eigenvector Method')
221
plt.legend()
222
plt.grid(True)
223
plt.xlim(0.2, 0.3)
224

225
plt.tight_layout()
226
plt.show()
227
```
228

229
### Model Order Selection
230

231
```python
232
import spectrum
233
import numpy as np
234
import matplotlib.pyplot as plt
235

236
# Generate test signal with known number of components
237
N = 100
238
n = np.arange(N)
239
frequencies = [0.1, 0.15, 0.3]  # Three sinusoids
240
amplitudes = [1.0, 0.8, 0.6]
241

242
signal = sum(A * np.exp(1j * 2 * np.pi * f * n) 
243
            for A, f in zip(amplitudes, frequencies))
244
noise = 0.1 * (np.random.randn(N) + 1j * np.random.randn(N))
245
noisy_signal = signal + noise
246

247
# Test different model orders
248
orders = range(5, 20)
249
aic_values = []
250
mdl_values = []
251

252
for order in orders:
253
    # Compute correlation matrix  
254
    R = spectrum.corrmtx(noisy_signal, order-1, method='autocorrelation')
255
    eigenvals = np.linalg.eigvals(R)
256
    eigenvals = np.sort(eigenvals)[::-1]  # Descending order
257
    
258
    # Compute criteria
259
    aic_vals = spectrum.aic_eigen(eigenvals, N)
260
    mdl_vals = spectrum.mdl_eigen(eigenvals, N)
261
    
262
    aic_values.append(np.min(aic_vals))
263
    mdl_values.append(np.min(mdl_vals))
264

265
# Plot model order selection
266
plt.figure(figsize=(10, 6))
267
plt.subplot(1, 2, 1)
268
plt.plot(orders, aic_values, 'bo-', label='AIC')
269
plt.xlabel('Model Order')
270
plt.ylabel('AIC Value')
271
plt.title('AIC Model Order Selection')
272
plt.grid(True)
273
optimal_aic = orders[np.argmin(aic_values)]
274
plt.axvline(optimal_aic, color='r', linestyle='--', 
275
           label=f'Optimal: {optimal_aic}')
276
plt.legend()
277

278
plt.subplot(1, 2, 2)  
279
plt.plot(orders, mdl_values, 'ro-', label='MDL')
280
plt.xlabel('Model Order')
281
plt.ylabel('MDL Value')
282
plt.title('MDL Model Order Selection')
283
plt.grid(True)
284
optimal_mdl = orders[np.argmin(mdl_values)]
285
plt.axvline(optimal_mdl, color='b', linestyle='--', 
286
           label=f'Optimal: {optimal_mdl}')
287
plt.legend()
288

289
plt.tight_layout()
290
plt.show()
291

292
print(f"True number of signals: {len(frequencies)}")
293
print(f"AIC optimal order: {optimal_aic}")
294
print(f"MDL optimal order: {optimal_mdl}")
295
```
296

297
### Minimum Variance Spectral Analysis
298

299
```python
300
import spectrum
301
import numpy as np
302
import matplotlib.pyplot as plt
303

304
# Generate signal with multiple frequency components
305
N = 128
306
t = np.arange(N)
307
fs = 1.0  # Sampling frequency
308

309
# Multiple sinusoids with different amplitudes
310
signal = (2.0 * np.sin(2*np.pi*0.1*t) +   # Strong component
311
         1.0 * np.sin(2*np.pi*0.25*t) +    # Medium component  
312
         0.5 * np.sin(2*np.pi*0.4*t))      # Weak component
313
noise = 0.2 * np.random.randn(N)
314
noisy_signal = signal + noise
315

316
# Compare different eigenvalue methods
317
orders = [10, 15, 20]
318
methods = ['MUSIC', 'Eigenvector', 'Minimum Variance']
319

320
plt.figure(figsize=(15, 10))
321

322
for i, order in enumerate(orders):
323
    plt.subplot(len(orders), 3, 3*i + 1)
324
    # MUSIC
325
    p_music = spectrum.pmusic(noisy_signal, IP=order, NSIG=3)
326
    freqs = np.linspace(0, fs/2, len(p_music.psd))
327
    plt.semilogy(freqs, p_music.psd)
328
    plt.title(f'MUSIC (Order={order})')
329
    plt.ylabel('Pseudospectrum')
330
    if i == len(orders)-1:
331
        plt.xlabel('Frequency')
332
    plt.grid(True)
333
    
334
    plt.subplot(len(orders), 3, 3*i + 2)
335
    # Eigenvector  
336
    p_ev = spectrum.pev(noisy_signal, IP=order, NSIG=3)
337
    plt.semilogy(freqs, p_ev.psd)
338
    plt.title(f'Eigenvector (Order={order})')
339
    if i == len(orders)-1:
340
        plt.xlabel('Frequency')
341
    plt.grid(True)
342
    
343
    plt.subplot(len(orders), 3, 3*i + 3)
344
    # Minimum Variance
345
    p_minvar = spectrum.pminvar(noisy_signal, order=order)
346
    plt.semilogy(freqs, p_minvar.psd)
347
    plt.title(f'Min Variance (Order={order})')
348
    if i == len(orders)-1:
349
        plt.xlabel('Frequency')
350
    plt.grid(True)
351

352
plt.tight_layout()
353
plt.show()
354
```