Отчет по лабораторной работе выполняется в виде связного (читаемого) текста в файле формата Microsoft

Лабораторная работа № 9

Исследование разомкнутой линейной системы при случайных возмущениях

• 
  освоение методов анализа одномерной линейной непрерывной системы при случайных возмущениях с помощью среды Matlab
  

• 
  научиться вычислять среднеквадратическое отклонение и дисперсию ан выходе линейной системы, возбуждаемой единичным белым шумом
  
• 
  научиться моделировать случайные процессы, используя в качестве источника сигнала белый шум (с ограниченной полосой)
  
• 
  научиться оценивать СКВО и дисперсию случайного процесса, полученного при моделировании
  
• 
  научиться вычислять автокорреляционную функцию случайного процесса
  
• 
  научиться вычислять спектральную плотность случайного процесса по известной корреляционной функции
  
• 
  научиться использовать быстрое преобразование Фурье (БПФ) для оценки спектральной плотности случайного процесса
  
• 
  научиться использовать спектральные окна для сглаживания оценки спектральной плотности случайного процесса
  

• 
  название предмета, номер и название лабораторной работы
  
• 
  фамилию и инициалы авторов, номер группы
  
• 
  фамилию и инициалы преподавателя
  
• 
  номер варианта
  
• 
  краткое описание исследуемой системы
  
• 
  результаты выполнения всех пунктов инструкции, которые выделены серым фоном (см. ниже): результаты вычислений, графики, ответы на вопросы.
  

Этап выполнения задания	Команды и иллюстрации
Очистите рабочее пространство Matlab (память).	clear all
Очистите окно Matlab.	clc
Введите передаточную функцию .	F = tf(1, [1 1])
Используя функцию norm, подсчитайте среднеквадратическое значение выхода этой системы при единичном белом шуме на входе.	norm ( F )
Подсчитайте дисперсию выхода системы при единичном белом шуме на входе.	norm ( F )^2
Найдите полосу пропускания этой системы (в рад/с).	bw = bandwidth ( F )
Найдите рекомендуемый максимальный интервал корреляции для моделирования по формуле	tau = 2*pi/100/bw
Запустите Simulink и создайте новую модель. Установите время моделирования 100 с (меню Simulation – Simulation Parameters – Stop Time).	на панели инструментов, в окне Simulink
Добавьте в модель блоки Band-Limited White Noise (белый шум с ограниченной полосой, группа Sources) и Scope (осциллограф, группа Sinks). Установите для белого шума параметр Noise Power (мощность) равный 1. Запустите модель и посмотрите, что представляет собой этот сигнал.
Так же, как и в предыдущей работе, подключите блоки Auto Correlator (автокорреляционная функция) и Power Spectral Density (спектральная плотность) из группы Simulink Extras – Additional Sinks). Посмотрите свойства этого сигнала.
Добавьте в схему звено с передаточной функцией так, как показано на схеме.
В параметрах блока Band-Limited White Noise уменьшите время корреляции (Sample Time) до значения, рассчитанного в п. 7. Для этого можно ввести в нужное поле имя переменной tau.
Откройте окно осциллографа, щелкните по кнопке и настройте параметры так, как показано на рисунке. На вкладке General в списке Sampling выберите вариант Sampling time (установить интервал вручную), а в соседнем поле введите имя переменной tau. На вкладке Data history нужно сбросить флажок Limit data points to last, установить флажок Save data to workspace, ввести имя переменной out и выбрать формат данных Array. Запустите моделирование.
Перейдите в окно Matlab, найдите среднеквадратическое отклонение (СКВО) и дисперсию сигнала на выходе звена. Сравните их со значениями, полученными в п. 4 и 5 по теоретическим формулам. Вычислите относительную ошибку при определении СКВО с помощью моделирования.	t = out(:,1); y = out(:,2); std ( y ) var ( y )
В окне блока Auto Correlator посмотрите, как выглядит корреляционная функция процесса, определенная по результатам моделирования.
Вычислите автокорреляционную функцию на выходе, используя функции Matlab, и постройте ее график. Сравните его с теоретической корреляционной функцией . Удобно создать новый m-файл и записать в него такие команды (без номеров строк): 1 R = xcorr(y)/length(y); 2 Rplus = R(floor(length(R)/2):end); 3 M = 200; 4 t = t(1:M); Rplus = Rplus(1:M); 5 R_teor = 0.5*exp(-abs(t)); 6 figure(1); 7 plot(t, Rplus, t, R_teor) 8 xlim([0 max(t)]); Комментарий: 1 – вычисляем экспериментальную корреляционную функцию, используя стандартную функцию xcorr из пакета Signal Processing 2 – выделяем корреляционную функцию для положительных 3 – число точек для оценки корреляционной функции (меньше, чем длина сигнала) 4 – «обрезаем» массив отсчетов времени и корреляционную функцию 5 – для тех же значений времени находим теоретическую корреляционную функцию 6-7 – построение обоих графиков на одном поле 8 – установить точные границы по оси абсцисс
Подключите блок Averaging Power Spectral Density (усредненная спектральная плотность) из группы Simulink Extras – Additional Sinks). Выполните моделирование еще раз и сравните спектры, полученные с помощью двух разных блоков. Сделайте выводы.
Постройте спектральную плотность сигнала для частот от 0 до 5 рад/с. По полученным данным. Сравните ее с теоретической спектральной плотностью . 1 T = t(2) - t(1); 2 w = 0:0.02:5; 3 Sw = []; Sw_teor =[]; 4 for i=1:length(w) 5 Sw(i) = sum(Rplus .* cos(w(i)*t)); 6 Sw_teor(i) = 1 / (w(i)*w(i) + 1); 7 end; 8 Sw = 2*T*Sw; 9 figure(2); 10 plot ( w, Sw, w, Sw_teor ); Комментарий: 1 – находим интервал дискретизации (он должен быть равен tau) 2 – задаем сетку частот, от 0 до 5 рад/с с шагом 0,02 рад/с 3 – опустошаем массивы 4-7 – цикл по всем выбранным частотам 5 – находим спектр как преобразование Фурье корреляционной функции 6 – теоретический спектр 8 – умножаем на 2T 9 – строим графики обоих спектров
Используя формулу , с помощью функции trapz (численное интегрирование методом трапеций), оцените дисперсию по экспериментальному и теоретическому спектрам. Объясните результаты, сравнив их со значением дисперсии, полученным в п. 5.	trapz(w,Sw)/pi trapz(w,Sw_teor)/pi
Постройте сглаженную оценку корреляционной функции с помощью окна Хэмминга. Для этого нужно добавить (в нужное место скрипта) команды hamm = 0.54 + 0.46*cos(pi*t/max(t)); Rhamm = Rplus .* hamm; и при построении корреляционных функций вывести третью линию plot(t, Rplus, t, R_teor, t, Rhamm) Перенесите в отчет полученный график.
Постройте оценку спектральной плотности, используя сглаженную корреляционную функцию. На графике должны быть три спектра (теоретический, оценка без сглаживания и оценка со сглаживанием). Скопируйте график в отчет.
Добавьте в скрипт команды для оценки спектральной плотности мощности с помощью быстрого преобразования Фурье (БПФ) и выполните его. 1 N = 2*pi/0.5/T; 2 N = 2^nextpow2(N); 3 Fw = T * fft(y, N); 4 Sw_fft = Fw .* conj(Fw) / N / T; 5 Sw_fft = Sw_fft(1:N/2+1); 6 w1 = 2*pi*[0:N/2] / N / T; 7 plot ( w, Sw_teor, w1, Sw_fft ); 8 xlim([0 max(w)]); Комментарий: 1 – считаем число точек для БПФ, чтобы шаг по частоте был равен 0,5 рад/с 2 – определяем ближайшую большую степень двойки 3 – выполняем БПФ 4 – считаем оценку спектра 5 – берем первую половину спектра до частоты Найквиста 6 – сетка угловых частот для построения графика 7 – строим теоретический спектр и оценку с помощью ДПФ 8 – устанавливаем пределы по оси абсцисс Сравните полученный результат с теоретической кривой. Сделайте выводы.
Повторите построение спектральной плотности, используя окно Хэмминга с масштабированием. Для этого добавьте в скрипт следующие команды: 1 scale = 1/sqrt(0.54^2 + 0.46^2/2); 2 hamm = hamming(N) * scale; 3 yHamm = y(1:N) .* hamm; 4 Fw = T * fft(yHamm, N); 5 Sw_fftHamm = Fw .* conj(Fw) / N / T; 6 Sw_fftHamm = Sw_fftHamm(1:N/2+1); 7 plot ( w, Sw_teor, w1, Sw_fft, w1, Sw_fftHamm ); 8 xlim([0 max(w)]); Комментарий: 1 – находим масштабирующий коэффициент для окна Хэмминга 2 – строим окно Хэмминга с масштабированием 3 – применяем окно к первым отсчетам сигнала 4 – выполняем БПФ 5 – считаем оценку спектра 6 – берем первую половину спектра до частоты Найквиста 7 – строим теоретический спектр и оценки 8 – устанавливаем пределы по оси абсцисс Запустите скрипт. Скопируйте полученный график в отчет. Сделайте выводы.