\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}

