Спектральна обробка сигналів: процедура та приклади застосування

ЛАБОРАТОРНА РОБОТА №4

Спектральна обробка сигналів: процедура та приклади застосування

Мета роботи:

1) практично ознайомитися з реалізацією процедури вимірювання спектру стаціонарного випадкового процесу в середовищі Matlab;

2)  змоделювати декілька прикладів застосування спектрального аналізу випадкових процесів, а саме:

а) виявлення періодичного сигналу, що маскується шумом;

б) вимірювання частот адитивної суміші кількох гармонічних коливань;

в) вимірювання частоти основного тону голосового сигналу.

Робоче завдання

1. Змоделювати задачу виявления періодичного сигналу, що маскується шумом, за умов: а) є апріорна інформація, що СВП являє собою адитивну суміш гармонічного процесу  з невідомими амплітудою , частотою  (значення якої знаходиться в межах 80-220 Гц), випадковою фазою, рівномірно розподіленою на інтервалі , та гаусівского білого шуму в смузі частот 0 – 5 кГц; б) відношення сигнал-шум цієї суміші дорівнює .
   • Змоделювати в середовищі Matlab адитивну суміш із заданими параметрами;
   • Розрахувати об’єм експериментальної вибірки відліків згенерованої суміші, що необхідна для забезпечення відношення сигнал-шум  на виході цифрового спектроаналізатора, що обчислює періодограму;
   • Обчислити та побудувати графіки різних оцінок спектру СВП : а) перідограма; б) оцінка Бартлета; в) оцінка Уелча із 50% перекриттям сегментів; г) авторегресійний метод Юла—Уолкера
   • Обчислити та побудувати графіки тих самих оцінок спектру СВП за умов збільшення  в 100 разів, при довжині сегменту , що дорівнює “старому” значенню
2. Здійснити спектральний аналіз сумиіші двох гармонічних коливань (без фонового шуму) за умов: а) частота першого коливання , частота другого коливання , де ; б) потужності обох гармонік однакові; в) =1024.
3. Здійснити натурний експеримент з вимірювання частоти основного тону голосового сигналу із застосуванням спектрального аналізу.
   • Використовуючи попередні записи слова та його фрагменту, ввести їх в робочий простір програми sptool;
   • здійснити спектральний аналіз із наступним вимірюванням частоти основного тону: а) перідограма; б) оцінка Бартлета; в) оцінка Уелча із 50% перекриттям сегментів; г) авторегресійний метод Юла—Уолкера

Таблиця 1. Варіанти значень числових параметрів

Оформити звіт за результатами експериментальних досліджень та аналітичних розрахунків за пп.1-3 даного робочого завдання (припустимо оформляти один звіт на бригаду).

1. Введення числових значень параметрів

Rvh = -13; f0 = 140; B = 5000;      % 4-й варіант

2. Моделювання задачі виявлення періодичного сигналу, що маскується шумом

Для адитивної суміші

(1)

сигналу  та шуму  вхідне відношення сигнал-шум:

,

або у децибелах:

.

Якщо приймемо , тоді параметр  дорівнює:

.

Відповідні команди Matlab:

%== обчисл.амплітуди А сигналу Y(t) ========

%

Dksi = 1;              % дисперс.шуму

A = 10^((10*log10(2)+Rvh)/20)       % амплітуда сигналу (відобр.в командному вікні)

= 0.3166

В лабораторній роботі №3 було побудовано графік фрагменту процесу , на якому сигнал не спостерігається. Крім того, в тій же роботі було показано, що кореляційний аналіз процесу  дозволяє виявити періодичний сигнал  на тлі шуму .

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

Оскільки складові частини процесу  статистично незалежні, враховуючи співвідношення Вінера-Хінчіна

,

одержуємо

,                                                                   (2)

де

,               (3)

– верхня гранична частота шуму .

Графік спектру суміші (2) побудувати неможливо, оскільки спектр необмеженого в часі гармонічного сигналу складається з двох дельта-функцій, що приймають безкінечно великі значення при  (рис.1).

На практиці завжди маємо відрізок, довжиною , реалізації випадкового процесу. Тому замість спектру потужності на практиці оперуємо із оцінкою спектру:

,

.

Можна показати, що математичне чекання відповідної оцінки відрізку гармонічного процесу має вигляд:

(4)

Тепер маємо ситуацію, максимально наближену до реальності, що дозволяє обчислити відношення сигнал-шум на виході спектроаналізатору. Якщо за сигнал на виході спектроаналізатору прийняти висоту піку математичного чекання (4), а за шум – математичне чекання рівня спектру шуму (насправді за рівень шуму приймають не рівень спектру шуму, а середньоквадритчну похибку оцінки спектра шуму – але можна показати, що для оцінки у вигляді періодограми ці два різні визначення відношення сигнал-шум співпадають), маємо:

,                                                                            (5)

тобто виграш у відношенні сигнал-шум за рахунок спектрального аналізу сягає величини

.                                                                               (6)

Звідси

,

.

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

Побудуємо графіки оцінки спектру суміші  та її математичного чекання. Зробимо це за два кроки.

Крок 1. Формування адитивної суміші гармонічного процесу та білого шуму

Rvyh = 10;

TB = 10^(0.1*(Rvyh-Rvh));

N = ceil(2*TB)     % кількість відліків процесу

N = 400

%== сигнал плюс шум ===============

%

fi = rand(1)*2*pi;             % випадкова початкова фаза

f0B = f0/B;                  % відносна частота сигналу

i = 1: N;

S = A*cos(pi*f0B*i+fi);         % N відліків сигналу

ksi = randn(1, N);                % N відліків шуму

Sksi = S + ksi;                       % N відліків суміші

subplot(3,1,1); plot(i,S);           % графік сигналу

title(‘Сигнал’);

subplot(3,1,2); plot(i,ksi);          % графік шуму

title(‘Шум’);

subplot(3,1,3); plot(i,Sksi,’r’);         % графік суміші

title(‘Сигнал плюс шум’);

Графіки процесів показано на рис.1.

Рис.1

Крок 2. Обчислення періодограми.

%=== періодограма =========

figure;

nfft = N; fs = 2*B; window = boxcar(N);

[Psksi,f1] = periodogram(Sksi,window,nfft,fs);

subplot(4,1,1); plot(f1,Psksi);

title(‘періодограма, лінійн.масштаб’);

subplot(4,1,3); plot(f1,10*log10(Psksi));

title(‘періодограма, логар.масштаб’);

% ==== матем.чекан. =========

f = 0:2*B/N:B;

Pksi = 1/B*rectpuls(f,2*B);

Ps = A^2*N/(4*B)*(sinc((f-f0)*N/2/B)).^2;

Psum = Pksi + Ps;

subplot(4,1,2); plot(f, Psum,’r’);

title(‘мат.чекан., лінійн.масштаб’);

subplot(4,1,4); plot(f, 10*log10(Psum),’r’);

title(‘мат.чекан., логар.масштаб’);

На рис.2 показані графіки оцінок спектрів та їх математичних чекань.

Рис.2

Розглянемо тепер такі різновиди оцінок спектру, що дають змогу зменшити дисперсію оцінки спектру білого шуму – це: а) оцінка Бартлета; в) оцінка Уелча із 50% перекриттям сегментів; г) авторегресійний метод Юла—Уолкера.

Вагові функції (вікна)

В даній лабораторній роботі для оцінки Уелча вам слід використати свій варіант вікна, що вказаний в таблиці 1.

Таблиця 1.

%=== оцінки спектру =====

% ==== періодограма =========

figure;

nfft = N; fs = 2*B; window = boxcar(N);

[Pmod,f1] = periodogram(Sksi,window,nfft,fs);

subplot(4,1,1); plot(f1,Pmod);

title(‘періодограма, лін.масшт.’);

%======= оцінка Бартлета =======

N1 = ceil(N/8);

nfft = N1; window = boxcar (N1); noverlap = 0;

[Pbart,f1] = pwelch(Sksi,window,noverlap,nfft,fs);

subplot(4,1,2); plot(f1,Pbart);

title(‘оц.Бартлета, лін..масшт.’);

%======= оцінка Уелча =======

N1 = ceil(N/8);

nfft = N1; Window = blackman (N1, ‘symmetric’); noverlap = ceil(N1/2);  % вікно Блекмана.

[Pwelc,f1] = pwelch(Sksi,window,noverlap,nfft,fs);

subplot(4,1,3); plot(f1,Pwelc);

title(‘оц. Уелча, лін..масшт.’);

%======= оцінка Юла-Уолкера =======

nfft = N;

p = 40;

[Pyul,f1] = pyulear(Sksi,p,nfft,fs);

subplot(4,1,4); plot(f1,Pyul);

title(‘оц. Юла-Уолкера, лін..масшт.’);

Результати спектрального оцінювання показані на рис.3.

Рис.3

Збільшимо  в 100 разів, а довжину сегменту  зробимо рівній “старому” значенню :

%== сигнал плюс шум ===============

%

fi = rand(1)*2*pi;             % випадкова початкова фаза

f0B = f0/B;                  % відносна частота сигналу

N=100*N;

i = 1:N;

S = A*cos(pi*f0B*i+fi);         % N відліків сигналу

ksi = randn(1, N);                % N відліків шуму

Sksi = S + ksi;                       % N відліків суміші

%=== оцінки спектру =====

% ==== періодограма =========

figure;

nfft = N; fs = 2*B; window = boxcar(N);

[Pmod,f1] = periodogram(Sksi,window,nfft,fs);

subplot(4,1,1); plot(f1,Pmod);

title(‘періодограма, лін.масшт.’);

%======= оцінка Бартлета =======

N1 = ceil(N/100);

nfft = N1; window = boxcar (N1); noverlap = 0;

[Pbart,f1] = pwelch(Sksi,window,noverlap,nfft,fs);

subplot(4,1,2); plot(f1,Pbart);

title(‘оц.Бартлета, лін..масшт.’);

%======= оцінка Уелча =======

N1 = ceil(N/100);

nfft = N1; Window = blackman (N1, ‘symmetric’); noverlap = ceil(N1/2);  % Блекмана.

[Pwelc,f1] = pwelch(Sksi,window,noverlap,nfft,fs);

subplot(4,1,3); plot(f1,Pwelc);

title(‘оц. Уелча, лін..масшт.’);

%======= оцінка Юла-Уолкера =======

nfft = N;

p = 40;

[Pyul,f1] = pyulear(Sksi,p,nfft,fs);

subplot(4,1,4); plot(f1,Pyul);

title(‘оц. Юла-Уолкера, лін..масшт.’);

Результати оцінювання наведені на рис.4.

Рис.4

3. Спектральний аналіз суміші двох гармонічних коливань

N = 1024;

f01 = f0;

f02 = f0 + 5*(2*B)/N;

k = 1:N;

Sksi = cos(2*pi*f01*k/(2*B))+ cos(2*pi*f02*k/(2*B));

% ==== періодограма =========

figure;

nfft = N; fs = 2*B; window = boxcar(N);

[Pmod,f1] = periodogram(Sksi,window,nfft,fs);

subplot(4,1,1); plot(f1(1:48),Pmod(1:48));

title(‘періодограма, лін.масшт.’);

%======= оцінка Бартлета =======

N1 = ceil(N/8);

nfft = N1; window = boxcar (N1); noverlap = 0;

[Pbart,f1] = pwelch(Sksi,window,noverlap,nfft,fs);

subplot(4,1,2); plot(f1(1:6),Pbart(1:6));

title(‘оц.Бартлета, лін..масшт.’);

%======= оцінка Уелча =======

N1 = ceil(N/8);

nfft = N1; Window = blackman (N1, ‘symmetric’); noverlap = ceil(N1/2);

[Pwelc,f1] = pwelch(Sksi,window,noverlap,nfft,fs);

subplot(4,1,3); plot(f1(1:6),Pwelc(1:6));

title(‘оц. Уелча, лін..масшт.’);

%======= оцінка Юла-Уолкера =======

nfft = N;

p = 40;

[Pyul,f1] = pyulear(Sksi,p,nfft,fs);

subplot(4,1,4); plot(f1(1:48),Pyul(1:48));

title(‘оц. Юла-Уолкера, лін..масшт.’);

Результати спектрального аналізу показані на рис.5.

Рис.5

4. Вимірювання частоти основного тону голосового сигналу із застосуванням спектрального методу

Активізуйте програму sptool та експортуйте в її робочий простір сигнал, що відповідає цілому слову та окремому голосному звуку із нього.

Обчисліть (в середовищі програми sptool) ті ж чотири різновиди оцінок спектру: періодограму, оцінку Бартлета, оцінку Уелча, оцінку Юла-Уолкера.

Параметр Nfft для оцінки у вигляді періодограми обчисліть за формулою:

,

де  і  – моменти початку та кінця фрагменту (=0); =11025 Гц.

Наприклад, для варіанту 4 було одержано =0.202 с, тому Nfft = 2227.

Приклад періодограми для варіанту 8 наведено на рис.6. На цьому рисунку бачимо два спектральних піки на частотах 138.6 Гц та 277 Гц – таким чином, оцінка частоти основного тону дорівнює 138.6 Гц

Рис.6

Оцінку Барлета для значення параметру   показано на рис.7. На цьому рисунку бачимо два спектральних піки на частотах 129 Гц та 258 Гц – таким чином, оцінка частоти основного тону дорівнює 129 Гц.

Рис.7

Оцінку Уелча для значення параметру , з вікном Кайзера та 50% перекриттям  показано на рис.8. На цьому рисунку бачимо два спектральних піки на частотах 129 Гц та 258 Гц – таким чином, оцінка частоти основного тону дорівнює 129 Гц.

Рис.8

Оцінку Юла-Уолкера для значення параметру Nfft = 2227 та порядку = 50 показано на рис.9 (для значень порядку = 80). На цьому рисунку бачимо два спектральних піки на частотах 138.6 Гц та 287 Гц – таким чином, оцінка частоти основного тону дорівнює 138.6 Гц.

Рис.9

Контрольні запитання: прокоментуйте одержані результати.

Додаток 2. Стислі теоретичні відомості

1. Непараметричні оцінки спектру

Періодограмою (periodogram) називається оцінка спектральної густини потужності:

Якщо при обчисленнях спектру викристовується вагова функція (вікно) , одержана оцінка спектру називається  модифікованою періодограмою (modified periodogram):

,

Можна показати, що періодограма не є слушною оцінкою спектру, оскільки дисперсія такої оцінки співставима їз квадратом її математичного чекання при довільному N:

Бартлетт (Bartlett) запропонував разділяти сигнал на сегменти, що не перекрываються, обчислювати для кожного сегменту періодограму і потім ці періодограми усереднювати (рис.10).

амплітуда

частота

час                                                              оцінка Бартлета

Рис.10

Уэлч (Welch) вніс в метод Бартлета два удосконалення: використання вагової функції й розбиття сигналу на сегменти, що перекриваються. Застосування вагової функції дозволяє ослабити розтікання спектру й зменшити зміщення оценки спектру ціною незначного погіршення разрішаючої здатності. Перекриття сегментів дозволяє зменшити дисперсію оцінки.

2. Обчислення непараметричних оцінок спектру в Matlab’і

В пакеті Signal Processing програми Matlab реалізовані три методи непараметричного оцінювання спектру — періодограма, метод Уелча й метод Томсона (MTM).

Обчислення періодограмми (в тому числі модифікованої) здійснюється за допомогою функції periodogram, синтаксис котрої має вигляд:

[Pxx,f] = periodogram(x,window,nfft,fs)

де вхідні дані – масив x відліків СВП, масив відліків вікна window, число точок БПФ nfft, частота дискретизації fs; вихідні дані – масиви значень перідограми Pxx та її аргумента f.

Метод Уелча є найбільш популярним періодограмним методом спектрального аналізу. В пакете Signal Processing він реалізується за допомогою функції pwelch:

[Pxx,f] = pwelch(x,window,noverlap,nfft,fs),

де noverlap – число відліків, що перекрываються.

3. Параметричні методи

Використання параметричних методів означає наявність математичної моделі випадкового процесу, що аналізується, а саме: СВП представляють як результат проходження білого шуму через фільтр (рис.11).

Білый шум                                                                               Забарвлений шум

Рис.11

Оскільки спектри СВП вхідного й вихідного процесів лінійної системи пов’язані співвідношенням:

,

для випадку білого шуму на вході, тобто , маемо

.

Таким чином, задачу вимірювання спектру потужності процесу  можна підмінити задачею вимірювання АЧХ деякого фільтру й потужності  білого шуму  на його вході.

Оскільки цифрові фільтри прийнято ділити на три типи, – нерекурсивні, рекурсивні й авторегресивні, – одержуєм три моделі процесу :

• ковзаючого середнього (КС);
• авторегресії ковзаючого середнього (АРКС);
• авторегресії (АР).

Спектральний аналіз зводиться в даному випадку до рішення оптимізаційної задачі, тобто до пошуку таких параметрів моделі, при яких вона найбільш близька до реального сигналу. В Matlab, в пакеті Signal Processing, реалізовані декілька авторегресійних оцінок спектру, а також оцінки, що базуються ще на двух методах, – на аналізі власних чисел й векторів кореляційної матриці сигналу: MUSIC (MUltiple SIgnal Classification) і EV (EigenVectors). В данній лабораторній роботі ми знайомимося тільки із авторегресійними оцінками спектру потужності ССП.

4. Реалізація обчислень в пакеті Matlab

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

Синтаксис функцій спектрального аналізу дуже схожий. Для авторегресійної оцінки Юла-Уолкера:

[Pxx,f] = pyulear(x,p,nfft,fs)

де p – порядок моделі. Інші позначення співпадають з такими для непе

5. Вимірювання частоти основного тону голосового сигналу із застосуванням спектрального методу

Активізуйте програму sptool та експортуйте в її робочий простір сигнал, що відповідає цілому слову та окремому голосному звуку із нього.

Обчисліть (в середовищі програми sptool) ті ж чотири різновиди оцінок спектру: періодограму, оцінку Бартлета, оцінку Уелча, оцінку Юла-Уолкера.

Параметр Nfft для оцінки у вигляді періодограми обчисліть за формулою:

,

де  і  – моменти початку та кінця фрагменту (=0); =11025 Гц.

Наприклад, для варіанту 4 було одержано =0.150 с, тому Nfft = 1654.

Періодограма для варіанту 4 наведено на рис.6. На цьому рисунку бачимо спектральний пік на частоті 219Гц – таким чином, оцінка частоти основного тону дорівнює 219 Гц

Рис.6

Оцінку Барлета показано на рис.7. На цьому рисунку бачимо спектральний пік на частоті 220 Гц  – таким чином, оцінка частоти основного тону дорівнює 220 Гц.

Рис.7

Оцінку Уелча для значення параметру , з вікном Блекмена та 50% перекриттям  показано на рис.8. На цьому рисунку бачимо спектральний пік на частотах 217 Гц – таким чином, оцінка частоти основного тону дорівнює 217 Гц.

Рис.8

Оцінку Юла-Уолкера для значення параметру Nfft = 1654 та порядку = 50 показано на рис.9. На цьому рисунку бачимо спектральний пік на частоті 219 Гц – таким чином, оцінка частоти основного тону дорівнює 219 Гц.

Рис.9

Висновки

В даній лабораторній роботі  ми практично ознайомилися з реалізацією процедури вимірювання спектру стаціонарного випадкового процесу в середовищі Matlab; впевнились у тому, що спектральний аналіз випадкового процесу допомагає вирішити задачу виявлення періодичного сигналу на тлі шуму. На практиці завжди маємо відрізок довжиною , тобто, обмеженого в часі, реалізації випадкового процесу. Тому замість спектру потужності на практиці оперуємо із оцінкою спектру.