0
# Non-parametric Methods
1

2
Classical spectral estimation techniques that do not assume an underlying parametric model. These methods are based on Fourier analysis and provide robust spectral estimates suitable for a wide range of signals and applications.
3

4
## Capabilities
5

6
### Periodogram Methods
7

8
Direct Fourier-based power spectral density estimation with various windowing and detrending options.
9

10
```python { .api }
11
def speriodogram(x, NFFT=None, detrend=True, sampling=1., scale_by_freq=True, window='hamming', axis=0):
12
    """
13
    Simple periodogram PSD estimation.
14
    
15
    Parameters:
16
    - x: array-like, input time series data
17
    - NFFT: int, FFT length (None uses len(x))
18
    - detrend: bool or str, detrending option (True, False, 'mean', 'linear')
19
    - sampling: float, sampling frequency
20
    - scale_by_freq: bool, scale PSD by frequency
21
    - window: str, window function name ('hamming', 'hanning', 'blackman', etc.)
22
    - axis: int, axis along which to compute PSD
23
    
24
    Returns:
25
    array: Power spectral density values
26
    """
27

28
class Periodogram(data, sampling=1., NFFT=None, window='hamming', detrend='mean', scale_by_freq=True):
29
    """
30
    Basic periodogram class for PSD estimation.
31
    
32
    Parameters:
33
    - data: array-like, input time series
34
    - sampling: float, sampling frequency
35
    - NFFT: int, FFT length for PSD computation
36
    - window: str, window function ('hanning', 'hamming', 'blackman', etc.)
37
    - detrend: str, detrending method ('mean', 'linear', None)
38
    - scale_by_freq: bool, scale PSD by frequency
39
    """
40
```
41

42
### Welch's Method
43

44
Averaged periodogram method that reduces variance by averaging multiple overlapping segments.
45

46
```python { .api }
47
def WelchPeriodogram(data, NFFT=None, sampling=1., **kargs):
48
    """
49
    Welch's averaged periodogram method.
50
    
51
    Parameters:
52
    - data: array-like, input time series data
53
    - NFFT: int, segment length for averaging
54
    - sampling: float, sampling frequency
55
    - **kargs: additional parameters (window, overlap, detrend)
56
    
57
    Returns:
58
    tuple: (frequencies array, PSD values array)
59
    """
60
```
61

62
### Daniell Periodogram
63

64
Periodogram with frequency domain smoothing using Daniell kernel for variance reduction.
65

66
```python { .api }
67
def DaniellPeriodogram(data, P, NFFT=None, detrend='mean', sampling=1., window_params={}):
68
    """
69
    Daniell's periodogram with smoothing.
70
    
71
    Parameters:
72
    - data: array-like, input time series data
73
    - P: int, smoothing parameter (Daniell kernel width)
74
    - NFFT: int, FFT length
75
    - detrend: str, detrending method
76
    - sampling: float, sampling frequency
77
    - window_params: dict, window function parameters
78
    
79
    Returns:
80
    tuple: (frequencies array, smoothed PSD array)
81
    """
82

83
class pdaniell(data, P, sampling=1., NFFT=None, detrend='mean', scale_by_freq=True):
84
    """
85
    Daniell periodogram class.
86
    
87
    Parameters:
88
    - data: array-like, input time series
89
    - P: int, Daniell smoothing parameter
90
    - sampling: float, sampling frequency
91
    - NFFT: int, FFT length for PSD computation
92
    - detrend: str, detrending method
93
    - scale_by_freq: bool, scale PSD by frequency
94
    """
95
```
96

97
### Correlogram Method
98

99
Blackman-Tukey method using windowed autocorrelation function for spectral estimation.
100

101
```python { .api }
102
103
class pcorrelogram(data, sampling=1., lag=-1, window='hamming', NFFT=None, scale_by_freq=True, detrend=None):
104
    """
105
    Correlogram PSD estimation class.
106
    
107
    Parameters:
108
    - data: array-like, input time series
109
    - sampling: float, sampling frequency
110
    - lag: int, maximum correlation lag
111
    - window: str, lag window function
112
    - NFFT: int, FFT length for PSD computation
113
    - scale_by_freq: bool, scale PSD by frequency
114
    - detrend: str, detrending method
115
    """
116
```
117

118
### Multi-Taper Method
119

120
Thomson's multi-taper method using discrete prolate spheroidal sequences for optimal spectral estimation.
121

122
```python { .api }
123
def pmtm(x, NW=None, k=None, NFFT=None, e=None, v=None, method="adapt", show=False):
124
    """
125
    Multi-taper method PSD estimation.
126
    
127
    Parameters:
128
    - x: array-like, input time series data
129
    - NW: float, time-bandwidth product (typically 2.5-4)
130
    - k: int, number of tapers (typically 2*NW - 1)
131
    - NFFT: int, FFT length
132
    - e: array-like, pre-computed tapers (DPSS sequences)
133
    - v: array-like, eigenvalues of tapers
134
    - method: str, averaging method ('adapt', 'unity', 'eigen')
135
    - show: bool, show intermediate plots
136
    
137
    Returns:
138
    tuple: (PSD array, frequencies array)
139
    """
140

141
def dpss(N, NW=None, k=None):
142
    """
143
    Discrete prolate spheroidal sequences (Slepian tapers).
144
    
145
    Parameters:
146
    - N: int, sequence length
147
    - NW: float, time-bandwidth product
148
    - k: int, number of sequences to compute
149
    
150
    Returns:
151
    tuple: (tapers array, eigenvalues array)
152
    """
153

154
class MultiTapering(Spectrum):
155
    """
156
    Multi-taper PSD estimation class.
157
    
158
    Parameters:
159
    - data: array-like, input time series
160
    - NW: float, time-bandwidth product
161
    - k: int, number of tapers
162
    - NFFT: int, FFT length
163
    - sampling: float, sampling frequency
164
    """
165
```
166

167
## Usage Examples
168

169
### Basic Periodogram Analysis
170

171
```python
172
import spectrum
173
import numpy as np
174
import matplotlib.pyplot as plt
175

176
# Generate test signal with multiple components
177
fs = 1000  # Sampling frequency
178
T = 1.0    # Duration
179
t = np.linspace(0, T, int(fs*T), False)
180

181
# Signal: two sinusoids + noise
182
f1, f2 = 50, 120  # Hz
183
signal = (np.sin(2*np.pi*f1*t) + 0.5*np.sin(2*np.pi*f2*t) + 
184
         0.2*np.random.randn(len(t)))
185

186
# Different periodogram methods
187
p_basic = spectrum.Periodogram(signal, sampling=fs)
188
p_welch = spectrum.WelchPeriodogram(signal, NFFT=256, sampling=fs)
189
p_daniell = spectrum.pdaniell(signal, P=5, sampling=fs)
190

191
# Plot comparison
192
plt.figure(figsize=(12, 10))
193

194
plt.subplot(3, 1, 1)
195
plt.plot(t[:200], signal[:200])
196
plt.title('Input Signal (first 200 samples)')
197
plt.xlabel('Time (s)')
198
plt.ylabel('Amplitude')
199
plt.grid(True)
200

201
plt.subplot(3, 1, 2)
202
freqs = np.linspace(0, fs/2, len(p_basic.psd))
203
plt.semilogy(freqs, p_basic.psd, label='Basic Periodogram')
204
plt.axvline(f1, color='r', linestyle='--', alpha=0.7, label=f'f1={f1}Hz')
205
plt.axvline(f2, color='r', linestyle='--', alpha=0.7, label=f'f2={f2}Hz')
206
plt.xlabel('Frequency (Hz)')
207
plt.ylabel('PSD')
208
plt.title('Basic Periodogram')
209
plt.legend()
210
plt.grid(True)
211
plt.xlim(0, 200)
212

213
plt.subplot(3, 1, 3)
214
freqs_welch, psd_welch = p_welch
215
plt.semilogy(freqs_welch, psd_welch, 'g-', label='Welch Method')
216
freqs_daniell = np.linspace(0, fs/2, len(p_daniell.psd))
217
plt.semilogy(freqs_daniell, p_daniell.psd, 'r-', label='Daniell Method')
218
plt.axvline(f1, color='k', linestyle='--', alpha=0.7)
219
plt.axvline(f2, color='k', linestyle='--', alpha=0.7)
220
plt.xlabel('Frequency (Hz)')
221
plt.ylabel('PSD')
222
plt.title('Welch vs Daniell Methods')
223
plt.legend()
224
plt.grid(True)
225
plt.xlim(0, 200)
226

227
plt.tight_layout()
228
plt.show()
229
```
230

231
### Multi-Taper Method Analysis
232

233
```python
234
import spectrum
235
import numpy as np
236
import matplotlib.pyplot as plt
237

238
# Generate signal with closely spaced frequencies
239
fs = 500
240
N = 500
241
t = np.arange(N) / fs
242

243
# Two closely spaced sinusoids + noise
244
f1, f2 = 100, 105  # 5 Hz apart
245
A1, A2 = 1.0, 0.8
246
noise_level = 0.3
247

248
signal = (A1*np.sin(2*np.pi*f1*t) + A2*np.sin(2*np.pi*f2*t) + 
249
         noise_level*np.random.randn(N))
250

251
# Multi-taper analysis with different parameters
252
NW_values = [2.5, 4.0, 6.0]
253
colors = ['blue', 'red', 'green']
254

255
plt.figure(figsize=(12, 10))
256

257
# Plot signal
258
plt.subplot(2, 2, 1)
259
plt.plot(t, signal)
260
plt.title(f'Signal: {f1}Hz + {f2}Hz + noise')
261
plt.xlabel('Time (s)')
262
plt.ylabel('Amplitude')
263
plt.grid(True)
264

265
# Compare with basic periodogram
266
plt.subplot(2, 2, 2)
267
p_period = spectrum.Periodogram(signal, sampling=fs)
268
freqs = np.linspace(0, fs/2, len(p_period.psd))
269
plt.semilogy(freqs, p_period.psd, 'k-', alpha=0.7, label='Periodogram')
270
plt.axvline(f1, color='r', linestyle=':', alpha=0.8)
271
plt.axvline(f2, color='r', linestyle=':', alpha=0.8)
272
plt.xlim(90, 115)
273
plt.xlabel('Frequency (Hz)')
274
plt.ylabel('PSD')
275
plt.title('Basic Periodogram')
276
plt.legend()
277
plt.grid(True)
278

279
# Multi-taper with different NW values
280
plt.subplot(2, 2, 3)
281
for i, NW in enumerate(NW_values):
282
    k = int(2*NW - 1)  # Number of tapers
283
    psd_mt, freqs_mt = spectrum.pmtm(signal, NW=NW, k=k, NFFT=1024)
284
    plt.semilogy(freqs_mt*fs, psd_mt, color=colors[i], 
285
                label=f'NW={NW}, k={k}', linewidth=2)
286

287
plt.axvline(f1, color='r', linestyle=':', alpha=0.8)
288
plt.axvline(f2, color='r', linestyle=':', alpha=0.8)
289
plt.xlim(90, 115)
290
plt.xlabel('Frequency (Hz)')
291
plt.ylabel('PSD')
292
plt.title('Multi-Taper Method (Different NW)')
293
plt.legend()
294
plt.grid(True)
295

296
# Taper visualization
297
plt.subplot(2, 2, 4)
298
NW = 4.0
299
k = int(2*NW - 1)
300
tapers, eigenvals = spectrum.dpss(N, NW=NW, k=k)
301
for i in range(min(4, k)):  # Show first 4 tapers
302
    plt.plot(tapers[:, i], label=f'Taper {i+1} (λ={eigenvals[i]:.3f})')
303
plt.xlabel('Sample')
304
plt.ylabel('Amplitude')  
305
plt.title(f'DPSS Tapers (N={N}, NW={NW})')
306
plt.legend()
307
plt.grid(True)
308

309
plt.tight_layout()
310
plt.show()
311

312
# Print taper properties
313
print(f"Multi-taper parameters: NW={NW}, k={k}")
314
print(f"Frequency resolution: {2*NW/N*fs:.2f} Hz")
315
print(f"Eigenvalue concentration: {np.sum(eigenvals > 0.9)}/{k} tapers > 0.9")
316
```
317

318
### Correlogram Method
319

320
```python
321
import spectrum
322
import numpy as np
323
import matplotlib.pyplot as plt
324

325
# Generate AR signal for correlogram analysis
326
N = 256
327
# AR(2) coefficients creating spectral peak around 0.2-0.3 normalized freq
328
ar_coeffs = [1, -1.6, 0.81]  # Poles near unit circle
329
noise = np.random.randn(N)
330

331
# Generate AR signal
332
from scipy import signal
333
ar_signal = signal.lfilter([1], ar_coeffs, noise)
334

335
# Different correlation lags for correlogram
336
lag_values = [20, 40, 80]
337
windows = ['rectangular', 'hamming', 'blackman']
338

339
plt.figure(figsize=(15, 10))
340

341
# Original signal
342
plt.subplot(2, 3, 1)
343
plt.plot(ar_signal)
344
plt.title('AR(2) Signal')
345
plt.xlabel('Sample')
346
plt.ylabel('Amplitude')
347
plt.grid(True)
348

349
# Autocorrelation function
350
plt.subplot(2, 3, 2)
351
max_lag = 50
352
autocorr = spectrum.CORRELATION(ar_signal, maxlags=max_lag, norm='unbiased')
353
lags = np.arange(-max_lag, max_lag+1)
354
plt.plot(lags, autocorr)
355
plt.title('Autocorrelation Function')
356
plt.xlabel('Lag')
357
plt.ylabel('Correlation')
358
plt.grid(True)
359

360
# Different lag lengths
361
plt.subplot(2, 3, 4)
362
for i, lag in enumerate(lag_values):
363
    p_corr = spectrum.pcorrelogram(ar_signal, lag=lag, window='hamming')
364
    freqs = np.linspace(0, 0.5, len(p_corr.psd))
365
    plt.semilogy(freqs, p_corr.psd, label=f'Lag={lag}', linewidth=2)
366
plt.xlabel('Normalized Frequency')
367
plt.ylabel('PSD')
368
plt.title('Correlogram: Different Lag Lengths')
369
plt.legend()
370
plt.grid(True)
371

372
# Different windows
373
plt.subplot(2, 3, 5)
374
lag = 40
375
for i, window in enumerate(windows):
376
    p_corr = spectrum.pcorrelogram(ar_signal, lag=lag, window=window)
377
    freqs = np.linspace(0, 0.5, len(p_corr.psd))
378
    plt.semilogy(freqs, p_corr.psd, label=f'{window.capitalize()}', linewidth=2)
379
plt.xlabel('Normalized Frequency')  
380
plt.ylabel('PSD')
381
plt.title(f'Correlogram: Different Windows (lag={lag})')
382
plt.legend()
383
plt.grid(True)
384

385
# Comparison with parametric method
386
plt.subplot(2, 3, 6)
387
p_corr = spectrum.pcorrelogram(ar_signal, lag=40, window='hamming')
388
p_burg = spectrum.pburg(ar_signal, order=10)
389
freqs = np.linspace(0, 0.5, len(p_corr.psd))
390
plt.semilogy(freqs, p_corr.psd, 'b-', label='Correlogram', linewidth=2)
391
plt.semilogy(freqs, p_burg.psd, 'r--', label='Burg AR(10)', linewidth=2)
392
plt.xlabel('Normalized Frequency')
393
plt.ylabel('PSD')
394
plt.title('Correlogram vs Parametric')
395
plt.legend()
396
plt.grid(True)
397

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

402
### Method Comparison for Different Signal Types
403

404
```python
405
import spectrum
406
import numpy as np
407
import matplotlib.pyplot as plt
408

409
def compare_methods(signal, title, fs=1.0):
410
    """Compare different non-parametric methods on a given signal."""
411
    
412
    # Method parameters
413
    methods = {
414
        'Periodogram': spectrum.Periodogram(signal, sampling=fs),
415
        'Correlogram': spectrum.pcorrelogram(signal, lag=len(signal)//4, 
416
                                           window='hamming'),
417
        'Multi-taper': None  # Will compute separately
418
    }
419
    
420
    # Multi-taper method
421
    NW = 3.0
422
    k = int(2*NW - 1)
423
    psd_mt, freqs_mt = spectrum.pmtm(signal, NW=NW, k=k, NFFT=512)
424
    
425
    plt.figure(figsize=(12, 8))
426
    
427
    # Plot signal
428
    plt.subplot(2, 2, 1)
429
    t = np.arange(len(signal)) / fs
430
    plt.plot(t, signal)
431
    plt.title(f'{title} - Time Domain')
432
    plt.xlabel('Time (s)')
433
    plt.ylabel('Amplitude')
434
    plt.grid(True)
435
    
436
    # Plot spectra
437
    plt.subplot(2, 2, 2)
438
    
439
    # Periodogram
440
    freqs = np.linspace(0, fs/2, len(methods['Periodogram'].psd))
441
    plt.semilogy(freqs, methods['Periodogram'].psd, 'b-', 
442
                label='Periodogram', linewidth=2)
443
    
444
    # Correlogram  
445
    freqs_corr = np.linspace(0, fs/2, len(methods['Correlogram'].psd))
446
    plt.semilogy(freqs_corr, methods['Correlogram'].psd, 'g--', 
447
                label='Correlogram', linewidth=2)
448
    
449
    # Multi-taper
450
    plt.semilogy(freqs_mt*fs, psd_mt, 'r:', 
451
                label='Multi-taper', linewidth=2)
452
    
453
    plt.xlabel('Frequency (Hz)')
454
    plt.ylabel('PSD')
455
    plt.title(f'{title} - PSD Comparison')
456
    plt.legend()
457
    plt.grid(True)
458
    
459
    plt.tight_layout()
460
    plt.show()
461

462
# Test different signal types
463
fs = 100
464

465
# 1. Narrow-band signal
466
N1 = 300
467
t1 = np.arange(N1) / fs
468
signal1 = np.sin(2*np.pi*10*t1) + 0.3*np.random.randn(N1)
469
compare_methods(signal1, "Narrow-band Signal (10 Hz)", fs)
470

471
# 2. Wide-band signal  
472
N2 = 400
473
signal2 = np.random.randn(N2)  # White noise
474
# Add some color with filtering
475
from scipy import signal as sp_signal
476
b, a = sp_signal.butter(4, 0.3, 'low')
477
signal2 = sp_signal.lfilter(b, a, signal2)
478
compare_methods(signal2, "Wide-band Colored Noise", fs)
479

480
# 3. Multi-component signal
481
N3 = 500
482
t3 = np.arange(N3) / fs  
483
signal3 = (np.sin(2*np.pi*5*t3) + 0.8*np.sin(2*np.pi*15*t3) + 
484
          0.6*np.sin(2*np.pi*25*t3) + 0.2*np.random.randn(N3))
485
compare_methods(signal3, "Multi-component Signal", fs)
486
```