Матриця ДПФ

У MATLAB розрахунок матриці прямого ДПФ реалізується за допомогою функції dftmtx. Синтаксис її виклику наступний:

А = dftmtx(n).

Тут n - розмірність ДПФ. Матриця зворотного ДПФ відрізняється від матриці прямого ДПФ комплексним сполученням і поділом на n, тому вона може бути отримана наступним чином: A_inv = conj(dftmtx(n))/n.

12.4. Блокова фільтрація в частотній області

Для реалізації блокової фільтрації за допомогою ШПФ в MATLAB використовується функція fftfilt, що має наступний синтаксис:

у = fftfilt(h, x, n)

Тут h - імпульсна характеристика фільтра, х - вхідний сигнал, n - розмірність ШПФ (точніше, розмірність ШПФ при завданні цього параметра визначається так: Nfft = 2^nextpow2(n), тобто n округляється вгору до ступеня двійки). Вихідний параметр у - результат фільтрації. Третій вхідний параметр n при виклику можна опускати, тоді розмірність ШПФ буде вибиратися автоматично, виходячи з максимальної ефективності обчислень (мінімального числа операцій).

12.5. Вікна

MATLAB містить (у пакеті Signal Processing) цілий ряд стандартних вагових функцій. Вони повертають вектори відліків, які можуть використовуватися в якості одного з параметрів різноманітних функцій непараметричного спектрального аналізу.

Всі розглянуті нижче функції приймають в якості параметра необхідну довжину вектора (n), яка повинна бути цілим позитивним числом, і повертають вектор-стовпець w. При n = 1 всі функції повертають значення 1. Графіки спектрів будуть будуватися функцією wvtool

Прямокутне вікно

Функція rectwin реалізує прямокутне вікно:

w = rectwin(n)

Повертаємий вектор заповнен одиницями: w = ones(n,1)

Візуалізатор вікон

Для одночасного перегляду віконних функцій у часовій і частотній областях призначений візуалізатор вікон, реалізований у вигляді функції wvtool (Window Visualization Tool). При її виклику необхідно вказати в списку вхідних параметрів вектор відліків вікна або кілька таких векторів:wvtool(w), wvtool (wl, w2, ...).

Як приклад викличемо функцію wvtool для перегляду прямокутного вікна 32 порядку:

>> wvtool(rectwin(32)) % На екрані з'явиться вікно для перегляду.

У лівій частині виводиться графік вікна в часовій області, а в правій - графік його амплітудного спектра. Під графіками відображається кількісна інформація про вікно:

Leakage factor — коефіцієнт витоку, який показує, яка частка загальної потужності вікна зосереджена в бічних пелюстках його спектра;

Relative sidelobe attenuation - рівень максимального з бічних пелюсток спектру щодо спектральної функції на нульовій частоті (тобто щодо головного пелюстка);

Mainlobe width (-3dB) — ширина головного пелюстка за рівнем -3 дБ. Тут мається на увазі повна ширина головного пелюстка (з урахуванням негативних частот), так що значення нормованої частоти, що на графіку відповідає рівню -3 дБ від максимуму, буде в два рази менше. Для наочності можна порівняти цю величину з нормованій відстанню між частотами сусідніх каналів ДПФ, що дорівнює 2/n (n - довжина вікна). У наведених прикладах це становить 2/32 = 0,0625.

Отже, прямокутне вікно забезпечує рівень бічних пелюсток -13,2 дБ, в бічних пелюстках спектру зосереджено 9,12% загальної потужності вікна, а нормована ширина головного пелюстка становить 0,055.

Трикутне вікно

Функція triang реалізує трикутне вікно: w = triangl(n).

>> wvtool(triang(32)) % На екрані з'явиться вікно для перегляду.

Косинусні вікна

Цей термін позначає вікна, які представляють собою суму декількох гармонічних складових, періоди яких укладаються на довжині вікна ціле число разів. До цієї категорії належать наступні вагові функції, реалізовані в пакеті Signal Processing:

Вікно Ханна - функція hann (за аналогією з вікном Хеммінга його часто помилково називають вікном Хеннінга; навіть відповідна функція в ранніх версіях пакета Signal Processing мала ім'я hanning):

w = hann(n, 'sflag').

>> wvtool(hann(32)) % На екрані з'явиться вікно для перегляду.

Вікно Хеммінга - функція: w = hamming(n, 'sflag').

>> wvtool(hamming(32)) % На екрані з'явиться вікно для перегляду.

Вікно Блекмена - функція: w = blackman(n, 'sflag').

>> wvtool(blackman(32)) % На екрані з'явиться вікно для перегляду.

Вікно Блекмена-Харріса - функція:w = blackmanharris(n).

>> wvtool(blackmanharris(32)) % На екрані з'явиться потрібне вікно.

Вікно Наттолла (версія вікна Блекмена-Харріса, запропонована Наттоллом) - функція: w = nuttallwin(n).

>> wvtool(nuttallwin(32)) % На екрані з'явиться вікно для перегляду.

Вікно з плоскою вершиною (flat top window) - функція: w = flattopwin(n, 'sflag').

>> wvtool(flattopwin(32)) % На екрані з'явиться вікно для перегляду.

Рядковий параметр 'sflag' дозволяє вибрати режим розрахунку вікна. При значенні 'symmetric', прийнятому за умовчанням, генерується симетричне вікно, для якого w(k) = w(n + 1 - k). При значенні 'periodic' створюється злегка несиметричне вікно, синусоїдальні компоненти якого будуть акуратно стикуватися при з'єднанні декількох примірників вікна.

Вікно Бомена

>> wvtool(bohmanwin(32)) % На екрані з'явиться вікно для перегляду.

Вікно Тьюкі

>> wvtool(tukeywin(32)) % На екрані з'явиться вікно для перегляду.

Вікно Кайзера

>> wvtool(kaiser(32,4)) % На екрані з'явиться вікно для β = 4.

>> wvtool(kaiser(32,9)) % На екрані з'явиться вікно для β = 9.

Вікно Чебишева

>> wvtool(chebwin(32,40)) % Для рівня бічних пелюсток – 40 дБ.

Гауссово вікно

>> wvtool(gausswin(32)) % На екрані з'явиться вікно для перегляду.

Вікно Парзена

>> wvtool(parzenwin(32)) % На екрані з'явиться вікно для перегляду.

12.6. Функції непараметричного спектрального аналізу

Розглянуто 3 функції, що здійснюють спектральний аналіз сигналу без використання додаткової інформації:

- specgram - обчислення миттєвого спектра сигналу,

- periodogram - обчислення спектральної щільності потужності однієї реалізації випадкового сигналу,

- pwelch - оцінка спектральної щільності потужності випадкового процесу методом усереднення модифікованих періодограмм.

Будуємо графік спектру щільності потужності дискретного випадкового процесу з експоненційної функцією кореляції.

function Example12_2

T=1; % період дискретизації

var_n=1; % дисперсія білого шуму

a1=0.9; % параметр експоненційної кореляції

f=linspace(0, 1/(2*T), 100); % вектор значень частот

P=2*var_n*T./abs(1-a1*exp(-i*2*pi*f*T)).^2;

hpsd = dspdata.psd(P, f); % Create a PSD data object.

plot(hpsd); % Plot the PSD.

load vcosig

specgram(vcosig, [], Fs) % Приклад використання функції specgram

12.7. Розрахунок періодограмми

Для обчислення періодограмми в MATLAB призначена функція periodogram. Синтаксис її виклику наступний:

[Рхх, f] = periodogram(x, window, Nfft, Fs, 'range')

Єдиним обов'язковим вхідним параметром є х - вектор відліків сигналу. Інші параметри мають значення за замовчуванням, які використовуються, якщо в якості параметра вказана порожня матриця [] або якщо деяка кількість параметрів (починаючи з останнього) опущені при виклику. Вектор window повинен містити коефіцієнти використовуваного вікна. За замовчуванням використовується прямокутне вікно. Параметр Nfft задає розмірність ШПФ, використовуваного для обчислення періодограмми. Параметр Fs - частота дискретизації в герцах. Рядковий параметр 'range' визначає частотний діапазон для повертаєтмого вектора Рхх. Можливі два значення:

□ 'twosided' - вектори Рхх і f мають довжину Nfft і відповідають повному діапазону частот 0...Fs. Цей варіант використовується за умовчанням, якщо х містить комплексні відліки;

□ 'onesided' - вектори Рхх і f мають довжину ceil((Nfft + 1)/2) і відповідають половинному діапазону частот 0...Fs/2. Цей варіант використовується за умовчанням в разі речового вектора х.

Якщо вихідні параметри при виклику не вказані, функція будує графік спектральної щільності за допомогою функції psdplot.

Як приклад оцінимо спектр щільності потужності експоненціально корельовано випадкового процесу.

function Example12_3

X0=randn(1, 1000); % формування випадкового сигналу

a=0.9;

X=filter(1, [1 -a], X0); % оцінка спектру щільності потужності

periodogram(X, [], [], 1)

На екрані періодограмма експоненціально корельовано випадкового процесу, отримана за допомогою функції periodogram

12.8. Функції авторегресійного спектрального аналізу

Синтаксис виклику функцій, що виробляють розрахунок спектру потужності на основі авторегресійної моделі формування сигналу, наведено нижче:

[Рхх, f] = рххх(х, p, Nfft, Fs, 'range')

Обов'язкові вхідні параметри: x - вектор відліків сигналу, р - порядок авторегресійної моделі. Решта вхідних параметрів мають значення за замовчуванням, які використовуються, якщо параметр представляє собою порожню матрицю [] або відсутній. Параметр Nfft задає число частотних точок для розрахунку спектра. Параметр Fs - частота дискретизації вхідного сигналу. Рядковий параметр 'range' визначає частотний діапазон для повертаємого вектора Рхх. Можливі два значення:

□ 'twosided' - вектори Рхх і f мають довжину Nfft і відповідають повному діапазону частот 0...Fs. Цей варіант використовується за умовчанням, якщо х містить комплексні відліки;

□ 'onesided' - вектори Рхх і f мають довжину ceil((Nfft + 1)/2) і відповідають половинному діапазону частот 0...Fs/2. Цей варіант використовується за умовчанням в разі речового вектора х.

В якості прикладу виконаємо спектральний аналіз експоненціально корельованого випадкового сигналу.

function Example12_4

X0=randn(1, 1000); % формування випадкового сигналу

a=0.9;

X=filter(1, [1 -a], X0); % оцінка спектру щільності потужності

pburg(X, 1, [], 1) % порядок моделі оптимальний

figure

pburg(X, 20, [], 1) % порядок моделі надмірний

На екрані результати авторегресійного спектрального аналізу при оптимальному і надмірному порядку моделі

Питання для самоперевірки

1. У чому відмінність прямого від зворотного ДПФ?

2. Навіщо потрібна функція fftshift?

3. Навіщо потрібні вікна і чому їх так багато?

4. Які функції непараметричного спектрального аналізу?

5. Як виконується розрахунок періодограмми?

6. Які функції авторегресійного спектрального аналізу?