\documentclass[a4paper,12pt]{extarticle} \usepackage{/home/eddy/ed} \title{Практикум \No2: статистика и вероятность} \author{}\date{}\nocolon \long\def\task#1{\noindent\leavevmode\refstepcounter{sect}\llap{\textbf{\thesect}\;}\indent\textit{#1}} \def\t#1{{\upshape\ttfamily #1}} \begin{document} \maketitle \section{Введение в Octave} Справка: \verb'help команда'. Случайные числа: \verb'rand(rows, cols)'~--- равномерное распределение, \verb'randn(rows, cols)'~--- нормальное, \verb'rande'~--- экспоненциальное, \verb'randg'~--- гамма-распределение, \verb'randp'~--- Пуассон. \subsection{Матрицы и векторы} Создание векторов: \verb'A=[1 2 3]' или же \verb'A=[1,2,3]'. Вектор-столбец: \verb|A=[1 2 3]'| или же \verb'A=[1;2;3]'. \verb'A=[начзн.:<шаг:>конзн.]'. Создание матриц: \verb'A=[1 5 2; 4 1 0]'. Конкатенация: \verb'A=[1:3]', \verb'B=[6:8]', \verb'C=[12:14]'; \verb'D=[A;B;C]'. Или: \verb|A=[1:3]'|, \verb|B=[6:8]'|, \verb|C=[12:14]'|; \verb'D=[A B C]'. Конкатенация возможна и для векторов и матриц: \verb|E=[(20:22)' (11:13)']|, \verb'[D E]'. \verb'X=[4 3 2]', \verb'Y=[9 8 7]', \verb'Z=[X Y]'. Размер в элементах: \verb'numel(A)'. Размер в байтах: \verb'sizeof(A)'. Размерность: \verb'size(A)'. Посмотреть значение переменной можно, указав ее имя или же команду \verb'disp(X)'. Специальные матрицы: \verb'eye(x)', \verb'ones(x)', \verb'zeros(x)', \verb'repmat(matrix, [ycount xcount])', \verb'magic(x)'. \textbf{Первая координата~--- Y!} Список переменных: \verb'whos'. Удалить переменные: \verb'clear'. \textbf{Индексы~--- с единицы!} Индексация: \verb'X=magic(5)', \verb'X(1:2,3:5)', \verb'X(2:end,3:end)', \verb'X(:,3)', \verb'X(3,:)', \verb'magic(4)(:,1)'. Матрицу~--- в столбец: \verb'magic(4)(:)'. Удаление элементов матриц и векторов: \verb'vec=[1:10]', \verb'idx=[3,7,5]', \verb'vec[idx]=[]'; \verb'X(:,[1,3])=[]'. \subsection{Матричные и векторные операции} Сложение, вычитание: \verb'A+B', \verb'A-B'. Транспонирование, сопряжение: \verb|A'|. Степень (для квадратных матриц): \verb'X^n'. \textbf{Отрицательная степень}: \verb'A^(-1)==inv(A)', но \verb'A^(-2)==inv(A)^2'! \textbf{Поэлементные операции}: \verb'A.(op)'. Умножение матриц: \verb'A*B', например, \begin{verbatim} A=[1 4 7; 8 5 2]; B=[4 5; 9 6; 3 1]; X=A*B ans = 61 36 83 72 \end{verbatim} \verb'A=[1 2; 3 4]', \verb'B=[50 60; 70 80]'. Матричное деление: \verb'A/B' (\verb'(A/B)*B==A'), эквивалентно $AB^{-1}$. <<Левое>> деление: \verb'A\B' (\verb'A*(A\B)==B'), эквивалентно $A^{-1}B$. \subsection{Графики} \verb'x=[0:0.1:2*pi]', \verb'plot(sin(x))'. Сохранить график: \verb'print -dpdf plot.pdf'. Гистограмма: \verb'hist(x,N)'. Например, \begin{verbatim} x=rand(1,10000); hist(x); hist(x, 5); x=randn(1,10000); hist(x,20); \end{verbatim} Графики двумерных функций. \verb'[X,Y]=meshgrid(-5:0.1:5,-6:0.1:6);', \verb'Z=X.^2-4*Y.^2;', \verb'surf(X,Y,Z)'. \verb'surfc(X,Y,Z)'~--- с изолиниями. \verb'mesh(X,Y,Z)'. Закрыть график: \verb'close all' или же \verb'close(fig)'. \verb'P=figure(N);'. Выбор между текущим графиком для отображения: \verb'figure(N)'. Например: \begin{verbatim} P1=figure(1); P2=figure(2); mesh(X,Y,Z); figure(1); surfc(X,Y,Z); close(P2); \end{verbatim} \verb'hold on/off'~--- дорисовать что-то на графике. Но интерфейс Octave к gnuplot не сравнится с самим gnuplot (пример)! Построим график плотности вероятности нормального распределения величины с $\mean{X}=-20$ и $\sigma_X=20$: \verb'x=[-70:30]; y=normpdf(x,-20,20);'. \subsection{Циклы, условия} Цикл \verb'for'. \verb'for x=1:10; printf("x=%d\n", x); endfor'. \begin{verbatim} X=[]; Y=[]; for x=1:3; A=dlmread(sprintf("for%02d", x)); > X=[X; A(:,1)]; Y=[Y; A(:,2)]; endfor plot(X, Y, 'o') \end{verbatim} Цикл \verb'while': \begin{verbatim} F=1; x=1; while F < 1e10; x+=1; F*=x; endwhile; printf("!%d=%g\n", x, F); !14=8.71783e+10 \end{verbatim} Условия: \begin{verbatim} x = input("Enter value: "); if(x < -5) printf("Less than -5\n") >elseif (x > 5) printf("More than 5\n") >else printf("Between -5 and 5\n") >endif \end{verbatim} \subsection{m-файлы} Функции и скрипты. Функции: \verb'checkX.m'. Скрипты: \verb'script.m'. Скрипт выполняется в глобальном пространстве! Пример (закомментировать \verb'x=' в \verb'script.m'): \begin{verbatim} clear x = [-2*pi:0.1:2*pi]; script \end{verbatim} Проверка переменной (и не только): \verb'exist(name, type)', скрипт \verb'script_chk' (запустить с определенной \verb'x' и после \verb'clear'). \subsection{Статистика} Среднее: \verb'mean(x)'. RMS: \verb'std(x)'. Медиана: \verb'median(x)'. Сумма: \verb'sum(x)'. Кумулятивная сумма: \verb'cumsum(x)'. Сортировка: \verb'sort(x)'. Минимум, максимум: \verb'min(x)', \verb'max(x)'. \section{Примеры выполнения заданий} \task{Сгенерировать синусоиду на участке~$[0,2\pi]$, добавить к ней гауссов белый шум с амплитудой~10\,дБ. Построить график.} Сгенерируем синусоидальный сигнал на участке~$[0,2\pi]$ командами $$\verb'x=[0:pi/50:2*pi]; y=sin(x);'$$ Теперь добавим к сигналу гауссов белый шум с амплитудой 10~дБ относительно амплитуды сигнала: $$\verb"y1=awgn(y,10,'measured'); plot(x, [y; y1])"$$ Третий параметр (measured) обязателен, т.к. без него процесс добавления шума будет несколько иным (мощность сигнала будет считаться равной 0~дБ), можете проверить на синусоиде с амплитудой~10. \task{Сгенерировать синусоиду с периодом $\pi/5$ на интервале $[0,20]$ с амплитудной модуляцией пилообразной функцией с периодом~10.} Для генерации синусоидального сигнала $y_0=A\sin(2\pi t/T)$ с амплитудной модуляцией по закону $y_1=f(t)$ необходимо перемножить эти две функции: $y=y_0\cdot y_1$. Промодулируем синусоиду с периодом $\pi/5$ пилообразным сигналом с периодом~10 на интервале $x\in[0,20]$. Для генерирования <<пилы>> используется функция \verb'sawtooth'. Если задать ей один аргумент (вектор $x$), период будет равен~$2\pi$, а сигнал будет изменяться в интервале $[-1,1]$. Чтобы задать смещение максимума, равное $a\cdot2\pi$, необходимо указать: \verb'y=sawtooth(x,a)'. Таким образом, чтобы получить <<пилу>> с интервалом сигнала в $[0,1]$ и периодом~10, необходимо дать команду \verb'y1=0.5+sawtooth(x*2*pi/10)/2;'. Следовательно, получить наш сигнал можно командой $$\verb'y=sin(x*10).*(0.5+sawtooth(x*pi/5)/2);'$$ (не забудьте про точку перед знаком умножения между функциями, иначе получите ошибку, т.к. Octave попытается перемножить два вектора--строки). \task{На интервале $[0,20]$ создать две синусоиды, сдвинутые друг относительно друга на 3~единицы. При помощи взаимно-корреляционной функции определить этот сдвиг.} \verb'x=[0:0.05:20];' \verb'y=sin(x);' \verb'y1=sin(x+3);' Попробуем определить, на сколько единиц сдвинут первый сигнал относительно второго. Для этого воспользуемся корреляционной функцией. Запишем \verb'Corr=xcorr(y,y1);'. Корреляционная функция в данном случае имеет вдвое большую ширину, чем исходная, т.к. она получается путем последовательного сдвига второй функции относительно первой. Поэтому построим график командой \verb'plot([-20:0.05:20],Corr)'. Воспользовавшись функцией увеличения можно увидеть, что ближайший к нулю максимум соответствует сдвигу одной функции относительно другой. В нашем случае сигнал был периодическим, поэтому при сдвигах на величины, превышающие половину периода, возникает ошибка, кратная полупериоду. Это необходимо учитывать в расчетах. Поменяем сдвиг на~0.2 и посмотрим, как изменится график. \task{Найдите сумму, разность, произведение и частное матриц $$A=\begin{pmatrix}4&12&8\\7&11&-1\\2&12&3\end{pmatrix},\qquad B=\begin{pmatrix}11&2&0\\0&3&0\\1&-1&12\end{pmatrix}.$$ Найдите определители исходных и получившихся матриц (команда \t{det(A)})} \begin{verbatim} a=[4 12 8; 7 11 -1; 2 12 3]; b=[11 2 0; 0 3 0; 1 -1 12]; sum=a+b dif=a-b prod=a*b div=a/b det(a) ... \end{verbatim} \section{Задания для самостоятельного выполнения} \begin{enumerate} \item Найдите сумму, разность, произведение и частное матриц $$A=\begin{pmatrix}1&2&3\\6&5&4\\9&8&7\end{pmatrix},\qquad B=\begin{pmatrix}9&8&7\\5&3&1\\0&2&6\end{pmatrix}.$$ Найдите определители исходных и получившихся матриц (команда \verb'det(A)'). \item Найдите значение почленного, матричного и скалярного произведений векторов $a=(2,5,7)$ и $b=(11,13,17)$. Скалярное произведение найдите двумя способами: путем перемножения векторов и при помощи функции \verb'dot(a,b)'. Найдите векторное произведение $a\times b$ при помощи функции \verb'cross(a,b)'. \item Постройте график нормального распределения на интервале $[0,100]$ с математическим ожиданием~50 и дисперсией~100. \item \label{noicy_AM} Получите сигнал с амплитудной модуляцией (из примера). Добавьте к нему гауссов белый шум с SNR 100, 50, 10 и~1~дБ. Постройте отдельно графики всех полученных сигналов. Можно ли сделать какой-либо вывод о виде сигнала при SNR=1? Как вы думаете, можно ли восстановить из него исходный сигнал? \item Для полученного сигнала найдите следующие характеристики: математическое ожидание (\verb'mean'), среднее квадратичное отклонение (\verb'std'), медиану (\verb'median') и моду (\verb'mode'). Найдите аналогичные величины для разности между зашумленным и оригинальным сигналом. Сравните полученные величины с теоретическими. \item Попытайтесь определить сдвиг двух синусоид (из примера) при зашумлении: \begin{itemize} \item только одной с уровнем сигнал/шум 1~дБ; \item обеих с уровнем SNR=1~дБ; \item одной с уровнем SNR=0.1~дБ; \item обеих с уровнем SNR=0.1~дБ. \end{itemize} Постройте один из сигналов с SNR=0.1~дБ. Можно ли определить его период? Можно ли определить период по автокорреляционной функции этого сигнала? \end{enumerate} \end{document}