Перед вами обучающий материал по основам цифровой обработки сигналов с использованием средств языка программирования Python. Предполагается, что читатель имеет базовые знания из области высшей математики, а также владеет языком Python и хотя бы поверхностно знает различные python-библиотеки - numpy/scipy, matplotlib и другие.
Для пользователей MATLAB / GNU Octave освоение материала с точки зрения программного кода не составит труда, поскольку основные функции и их атрибуты во многом идентичны и схожи с методами из python-библиотек.
Нередко перед исследователями встаёт вопрос анализа характеристик сигнала, процесса или явления.
Одним из инструментов по определению таких характеристик является спектральный анализ.
Методы спектрального анализа можно разделить на две большие группы:
Для применения параметрических методов требуется априорное знание о параметрах исследуемого объекта, в случае непараметрических методов - не требуется. Именно о последних и пойдет речь в данном семинаре.
Допущение: Будут рассматриваться стационарные в широком смысле (WSS - wide-sense stationary) случайные процессы.
Спектр мощности для случайного в широком смысле процесса $x(n)$ (где $n$ - номер временного отсчета) представляет из себя ничто иное, как преобразование Фурье автокорреляционной функции данного процесса [1, 393]:
$$ P_x(e^{j\omega}) = \sum_{-\infty}^{\infty} r_x(k)e^{jk\omega} \tag{9.1}$$где $\omega$ - это циклическая частота, а $k$ - номер частотного отсчета (индекс частоты).
На конечном интервале автокорреляционная функция примет вид:
$$ \hat{r}(k) = \frac{1}{N}\sum_{n=0}^{N-1}x(n+k)x^*(n) \tag{9.2}$$где $N$ - длина последовательности временных отсчетов.
Одним из самых простых и популярных представителей группы непараметрических методов является метод периодограмм, предложенный Артуром Шустером ещё в конце XIX века. Рассмотрим математическое описание данного метода.
Для начала, в качестве ограничения зададим условие, непозволяющее нашему сигналу $x(n)$ выходить за пределы интервала:
$$ x_N(n) = \begin{cases} x(n), &0 \leq n \leq N \\ 0, &\text{otherwise} \end{cases} \tag{9.3}$$Следовательно, $x_N(n)$ - это результат перемножения сигнала с оконной функцией прямоугольного вида: $$ x_N(n) = \omega_R(n)x(n) \tag{9.4}$$
Приняв во внимание формулу (3), переопределим формулу (2) через свертку: $$ \hat{r}(k) = \frac{1}{N}x_N(k)*x_N(-k) \tag{9.5}$$
Тогда преобразование Фурье автокорреляционной функции даст следующий результат: $$ \hat{P}_{Per}(e^{j\omega}) = \frac{1}{N}X_N(e^{j\omega})X_N^*(e^{j\omega}) = \frac{1}{N} \left| X_N(e^{j\omega}) \right|^2 \tag{9.6}$$
Следовательно, периодограмма пропорциональна квадрату амплитуды преобразованию Фурье дискретного времени (DTFT - discrete-time Fourier transform), а значит вполне может быть вычислена через алгоритм БПФ (FFT) на этапе программной реализации.
В python метод периодограмм реализован в рамках библиотеки scipy.signals в виде метода periodogram.
Возвратимся к формуле (4). Было отмечено, что периодограмма пропорциональна результату перемножения сигнала с оконной функцией $\omega(n)$, при условии, что данная функцию является прямоугольной. Однако, нужно отметить, что оконная функция может быть и других форм.
В методе periodogram библиотеки scipy форма окна определяется параметром window (полный список доступных форм можно посмотреть по данной ссылке).
Периодограмма с использованием непрямоугольного окна часто называется модифицированной:
$$ \hat{P}_{M}(e^{j\omega}) = \frac{1}{NU}\left| \sum_{n=-\infty}^{\infty} x(n)\omega(n)e^{-jn\omega}\right|^2 \tag{9.7}$$где $N$ - это длина окна, а $U$ - это константа [1, с. 410]: $$ U= \frac{1}{N}\sum_{n=0}^{N-1}|\omega(n)|^2 \tag{9.8}$$
которая показывает асимптотическую несмещенность (unbiased) модифицированной периодограммы.
Рассмотрим более последовательные оценки мощностных спектров. Предпосылкой является наблюдение, что, при увеличении длины последовательности $N$ до бесконечности, математическое ожидание периодограммы стремится к $P_x(e^{j\omega})$:
$$ \lim_{N \to \infty} E\{ \hat{P}_{Per}(e^{j\omega})\} = P_x(e^{j\omega}) \tag{9.9}$$Соответственно, если мы найдём последовательную оценку мат. ожидания $E\{ \hat{P}_{Per}(e^{j\omega})\}$, то оценка $P_x(e^{j\omega})$ так же будет последовательной. Для этого можно применить классическое усреднение по некоторой выбоке реализаций.
Конечно, в реальных системах более реалистичным является случай, когда вместо сбора достаточного количества реализаций процесса, мы исследуем сигнал достаточной длины, который разбивается на последовательности (sequences), а уже те, в свою очередь, используются для усреднения (рис.1).
Рис. 1. Разбиение $x(n)$ на неперекрывающиеся последовательности [1, c. 413].
В литературе такой подход называется методом Бартлетта [2, c. 332]:
$$ \hat{P}_B(e^{j\omega}) = \frac{1}{N} \sum_{i=0}^{K-1}\left| \sum_{n=0}^{L-1} x(n+iL)e^{-jn\omega} \right|^2 \tag{9.10}$$где $K$ - это количество неперекрывающихся (non-overlapping) последовательностей длинной $L$ каждая, а $N=KL$. Иначе говоря, перед нами формула усреднения периодограмм.
Однако, и это ещё не всё. В 1967 году Ф.Д. Уэлч предлагает метод, который позже будет носить его же имя [2, c.333][3]. Подразумеваются следующие отличия от метода Бартлетта:
разбиение сигнала в том числе на пересекающиеся (overlapping) последовательности;
применение не только прямоугольных оконных функций (модифицированных периодограмм).
В python метод Уэлча реализован в рамках библиотеки scipy.signals в виде метода welch.
Применение непрямоугольных окон позволяет достичь больших степеней свободы в вопросах перекрытия (overlapping) [4]. При этом, используя перекрывающиеся последовательности, мы повышаем общее количество фрагметов сигнала. Как следствие, периодограмма Уэлча будет менее осциллирующей (изрезанной), чем периодограмма Бартлетта [2, c.333].
Можно отметить следующую тенденции: метод Уэлча является общим случаем для методов Бартлетта и модифицированной периодограммы, а значит и для самого метода периодограмм. Именно это наблюдение легло в основу архитектуры методов в рамках реализации их на Python в рамках библиотеки scipy.
import numpy as np
from scipy import signal
import matplotlib.pyplot as plt
FONT_SMALL = 12
FONT_MEDIUM = 14
plt.rc('axes', titlesize=FONT_MEDIUM)
plt.rc('axes', labelsize=FONT_SMALL)
plt.rc('xtick', labelsize=FONT_SMALL)
plt.rc('ytick', labelsize=FONT_SMALL)
plt.rc('legend', fontsize=FONT_MEDIUM)
plt.rc('figure', titlesize=FONT_MEDIUM)
%matplotlib inline
В качестве примера рассмотрим гармонический сигнал следующего вида:
w_1 = 40 # frequency of the 1st component of the signal (Hz)
w_2 = 60 # frequency of the 2nd component of the signal (Hz)
a = 0.5 # magnitude of the 1st component of the signal
b = 1.0 # magnitude of the 2nd component of the signal
Зададим временной интервал:
Nsub = 10 # number of subsequences
t = np.array([i for i in range(1,301*Nsub)])/1000 # time samples (s)
fs = 1 / (t[1]-t[0]) # sampling frequency (Hz)
Nfft = int(3e3)
Моделируем сигнал:
x = a*np.cos(2*np.pi*w_1*t) + b*np.sin(2*np.pi*w_2*t) # considered signal
f = fs*np.array([i for i in range(int(len(x)))]) / len(x) # frequencies
# Plot results
plt.subplots(1, 1, figsize=(16, 4))
plt.plot(t[:200], x[:200], lw=1.5)
plt.ylabel('Signal magnitudes')
plt.xlabel('Time (s)')
plt.title('Signal')
plt.xlim([0, 0.2])
plt.grid(True)
plt.tick_params(axis ='x', rotation = 45)
plt.tight_layout()
Построим частотный спектр амплитуд ДПФ:
FFT = np.fft.fft(x, n=Nfft) # Fast Fourier Transform
amps = np.abs(FFT) / (len(FFT) / 2) # magnitudes of FFT
l = int(len(f)/4)
plt.subplots(1, 1, figsize=(16, 4))
plt.stem(f[:l], amps[:l], use_line_collection=True)
plt.ylabel('FFT magnitudes')
plt.xlabel('Frequencies (Hz)')
plt.title('Signal')
plt.grid(True)
plt.xticks(np.arange(0, f[l], 10))
plt.yticks(np.arange(0, max(amps)+0.1, .1))
plt.tick_params(axis ='x', rotation = 45)
plt.xlim([0, 240])
plt.tight_layout()
Моделируем аддитивный шум:
np.random.seed(42)
n = 2*np.random.randn(len(t)) # white Gaussian noise
FFT = np.fft.fft(n, n=Nfft)
amps = np.abs(FFT) / (len(FFT) / 2)
plt.subplots(1, 1, figsize=(16, 4))
plt.stem(f[:len(amps)], amps, use_line_collection=True)
plt.ylabel('FFT magnitudes')
plt.xlabel('Frequencies (Hz)')
plt.title('Noise')
plt.yticks(np.arange(0, max(amps), .1))
plt.xlim([0, 1000])
plt.grid(True)
plt.tight_layout()
Добавляем шум к исходному сигналу:
y = x + n # signal + noise
plt.subplots(1, 1, figsize=(16, 4))
plt.plot(t[:200], y[:200])
plt.ylabel('Noised signal magnitudes')
plt.xlabel('Time (s)')
plt.title('Signal + noise')
plt.xlim([0, 0.2])
plt.tick_params(axis ='x', rotation = 45)
plt.grid(True)
plt.tight_layout()
Попробуем проанализировать данный сигнал с помощью методов периодограмм и модифицированных периодограмм.
windows = ['hamming', None]
plt.subplots(1, 1, figsize=(16, 6))
for window in windows:
if window == None:
label = 'rectangular'
else:
label = window
f, Pxx_den = signal.periodogram(y, fs=fs, scaling='spectrum', nfft=Nfft, window=window)
plt.semilogy(f[1:], Pxx_den[1:], label=label)
plt.ylabel('Spectrum')
plt.xlabel('Frequencies (Hz)')
plt.title('Periodogram')
plt.xticks(np.arange(0, max(f), 10))
plt.xlim(0, 500)
plt.tick_params(axis ='x', rotation = 70)
plt.legend(loc='best')
plt.grid(True)
plt.tight_layout()
Две составляющие полезного сигнала, конечно, различимы, но при больших шумовых амплитудах методы могут и не справится.
Попробуем применить метод Бартлетта для того же сигнала:
windows = ['bartlett', 'hamming', 'blackman']
plt.subplots(1, 1, figsize=(16, 6), dpi=100)
for window in windows:
if window == 'bartlett':
label = 'rectangular'
else:
label = window
f, Pxx_den = signal.welch(y, fs=fs, nperseg = len(x)/Nsub, noverlap=0, scaling='spectrum', nfft=Nfft, window=window)
plt.semilogy(f, Pxx_den, label=label)
plt.legend()
plt.ylabel('Spectrum')
plt.xlabel('Frequencies (Hz)')
plt.title('Bartlett')
plt.xlim(0, 500)
plt.xticks(np.arange(0, max(f), 10))
plt.tick_params(axis ='x', rotation = 70)
plt.legend(loc='best')
plt.grid(True)
plt.tight_layout()
Отмечаем лучшую различимость составляющих полезного сигнала в случае разбиения на последовательности.