Допоміжні відомості для роботи з програмою Matlab
Введення числових значень параметрів
Для автоматизації виконання даної лабораторної роботи введіть в робочий простір програми Matlab числові значення параметрів відоповідно до Таблиці 1. Для цього скопіюйте строку із даними свого варіанту в командне вікно програми Matlab та натисніть клавішу Enter.
Rvh = -10; f0 = 80; B = 5000; % 1-й варіант
Rvh = -11 f0 = 100; B = 5000; % 2-й варіант
Rvh = -12; f0 = 120; B = 5000; % 3-й варіант
Rvh = -13; f0 = 140; B = 5000; % 4-й варіант
Rvh = -14; f0 = 160; B = 5000; % 5-й варіант
Rvh = -15; f0 = 180; B = 5000; % 6-й варіант
Rvh = -16; f0 = 200; B = 5000; % 7-й варіант
Rvh = -17; f0 = 220; B = 5000; % 8-й варіант
2. Моделюваня задачі виявления періодичного сигналу, що маскується шумом
Для адитивної суміші
(1)
сигналу
та шуму
вхідне відношення сигнал-шум:
,
або у децибелах:
.
Якщо
приймемо
,
тоді параметр
дорівнює:
.
Відповідні команди Matlab:
%== обчисл.амплітуди А сигналу Y(t) ========
%
Dksi = 1; % дисперс.шуму
A = 10^((10*log10(2)+Rvh)/20) % амплітуда сигналу (відобр.в командному вікні)
Для варіанту 8 одержуємо = 0.1998
В
лабораторній
роботі №3
було побудовано графік фрагменту процесу
,
на якому сигнал не спостерігається.
Крім того, в тій же роботі було показано,
що кореляційний аналіз процесу
дозволяє виявити періодичний сигнал
на тлі шуму
.
Покажемо тепер, что спектральний аналіз випадкового процесу допомагає вирішити задачу виявлення періодичного сигналу на тлі шуму.
Оскільки складові частини процесу статистично незалежні, враховуючи співвідношення Вінера-Хінчіна
,
одержуємо
,
(2)
де
,
(3)
-
верхня гранична частота шуму
.
Контрольне запитання: доведіть справедливість співвідношень (2) та (3), спираючись на відомі вам формули кореляційних функцій гармонічного сигналу та білого шуму, обмеженого за частотою.
Графік
спектру суміші (2) побудувати неможливо,
оскільки спектр необмеженого в часі
гармонічного сигналу складається з
двох дельта-функцій, що приймають
безкінечно великі значення при
(рис.1).
На
практиці завжди маємо відрізок, довжиною
,
реалізації випадкового процесу. Тому
замість спектру потужності на практиці
оперуємо із оцінкою
спектру:
,
.
Можна показати, що математичне чекання відповідної оценки відрізку гармонічного процесу має вигляд:
(4)
Тепер маємо ситуацію, максимально наближену до реальності, що дозволяє обчислити відношення сигнал-шум на виході спектроаналізатору. Якщо за сигнал на виході спектроаналізатору прийняти висоту піку математичного чекання (4), а за шум - математичне чекання рівня спектру шуму (насправді за рівень шуму приймають не рівень спектру шуму, а середньоквадритчну похибку оцінки спектра шуму – але можна показати, що для оцінки у вигляді періодограми ці два різні визначення відношення сигнал-шум співпадають), маємо:
,
(5)
тобто виграш у відношенні сигнал-шум за рахунок спектрального аналізу сягає величини
.
(6)
Звідси
,
.
Обрахувавши
та враховучи, що за теоремою Котельнікова
неперервний процес можна дискретизувати
з кроком
,
одержуємо корисну для планування
експерименту формулу, що дозволяє
обрахувати потрібну кількість відліків
процесу:
Побудуємо графіки оцінки спектру суміші та її математичного чекання. Зробимо це за два кроки.
Крок 1. Формування адитивної суміші гармонічного процесу та білого шуму
Rvyh = 10;
TB = 10^(0.1*(Rvyh-Rvh));
N = ceil(2*TB) % кількість відліків процесу
%== сигнал плюс шум ===============
%
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
Контрольні запитання: 1) чому спектр гармонічного сигналу видно недосить добре на тлі спектру шуму, хоча відношення сигнал шум дорівнює 10 дБ? 2) чому гармонічний сигнал гірше спостерігається на графіку періодограми в логаріфмічному масштабі?
Розглянемо тепер такі різновиди оцінок спектру, що дають змогу зменшити дисперсію оцінки спектру білого шуму – це: а) оцінка Бартлета; в) оцінка Уелча із 50% перекриттям сегментів; г) авторегресійний метод Юла—Уолкера.
Вагові функції (вікна)
В даній лабораторній роботі для оцінки Уелча вам слід використати свій варіант вікна, що вказаний в таблиці 1.
Таблиця 1.
№ варіанту |
Вікно |
Синтаксис |
1 |
Прямокутне |
Window = boxcar (N1) |
2 |
Трикутне |
Window = triang (N1) |
3 |
Бартлета |
Window = bartlett (N1) |
4 |
Блекмена |
Window = blackman (N1, ‘symmetric’) |
5 |
Чебишева (40db) |
Window = chebyshev (N1, 40) |
6 |
Хемінга |
Window = hamming (N1, ‘symmetric’) |
7 |
Хана (Хенінга) |
Window = hanning (N1, ‘symmetric’) |
8 |
Кайзера (beta = 6) |
Window = kaiser (N1, 6) |
%=== оцінки спектру =====
% ==== періодограма =========
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 = kaiser(N1,6); noverlap = ceil(N1/2); % вікно Кайзера з перекрит. 50%.
[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
Контрольні запитання: прокоментуйте рис.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 = kaiser(N1,6); noverlap = ceil(N1/2); % вікно Кайзера з перекрит. 50%.
[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
Контрольні запитання: прокоментуйте рис.4, у порівнянні із рис.3