Сведения о программировании в системе инженерных и научных расчетов MATLAB
MATLAB – среда для математических вычислений, разработки алгоритмов, визуализации и анализа данных.
Операторы
В среде MATLAB используются операторы +, -, *, /, =, ^ (возведение в степень).
Переменные
Переменные объявляются следующим образом:
N = 1024;
Z = 28.2;
Имя переменной может содержать только символы A...Z, a...z, 0...9 и _, и не должно начинаться с цифры или подчеркивания. Точка с запятой в конце строки блокирует вывод значения переменной в командное окно.
Комментарий
Для задания комментария используется символ %. Например:
% Комментарий
Функции
clear all - очищает рабочее пространство – удаляет все переменные раннее запущенной программы.
сlc -очищает командное окно.
Функция digits(n) указывает количество значимых цифр после запятой для расчетов. По умолчанию – 32.
Функции cos(X) и sin(Y) возвращают косинус и синус угла X, заданного в радианах.
Функция sqrt(X) возвращает квадратный корень числа X.
Функция fix(A) возвращает число A, округленное в меньшую сторону. Например, fix(5.7) возвращает число 5.
Функция rand(n) возвращает массив размера n x n, заполненный случайными числами в диапазоне от 0 до 1. Например, rand(1) возвращает одно случайное число в диапазоне от 0 до 1.
Циклы
Цикл по переменной i от 1 до N реализуется следующим образом:
for i = 1:N,
…
end
Пример цикла по переменной iи вложенного в него цикла по переменной j:
for i = 1:N1,
for j = 1:N2,
…
end
end
Массивы
Для выделения массива в памяти используется функция zeros.
B = zeros(n) создает массив B размера n x n, заполненный нулями.
B = zeros(m,n)или B = zeros([m n]) создает массив B размера m x n, заполненный нулями.
B = zeros(m,n,p,...)илиB = zeros([m n p ...])создает массив B размера m x n x p x…, заполненный нулями.
Пример выделения и заполнения массива:
B = zeros(1,180);
for i = 1:180,
B(i) = sin(i/180*pi);
end
Условие
Условие реализуется следующим образом:
if (A < B),
C = A + B;
else
C = A - B;
end
Если выполняется условие A < B, то C = A + B, иначе C = A - B.
Если требуется объединить два или более условий, то используются оператор && для логической операции «И» или оператор || для логической операции «ИЛИ». Например:
if ((A < B) &&(A > 5)),
Задание интервала
Код i = 1:N заполняет массив i значениями от 1 до N.
Графики
Для вывода графиков используются следующие функции:
figure -Выводит график в отдельном окне.
Функция subplot(mnp) разбивает окно на m x n графиков и выбирает p-й график как текущий.
Функция plot(X,Y,'Color') строит график массивов X по оси абсцисс и Y по оси ординат цветом линии Color.
| Значение параметра Color | Цвет |
| r | Красный |
| g | Зеленый |
| b | Синий |
| c | Голубой |
| m | Розовый |
| y | Желтый |
| k | Черный |
title(' ')выводит название графика.
xlabel(' ') выводит текст по оси абсцисс.
ylabel(' ')выводит текст по оси ординат.
Функция xlim([xmin xmax]) указывает диапазон значений вывода по оси абсцисс от xmin до xmax.
Функция ylim([ymin ymax]) указывает диапазон значений вывода по оси ординат от ymin до ymax.
Пример программы моделирования цифровой фильтрации зашумленного гармонического сигнала:
% Фильтрация зашумленного гармонического сигнала
clc; % очистить командное окно
clear all; % очистить рабочее пространство
digits(4); % точность расчетов 4 знака после запятой
%--------------------------------------------------------------------
N = 1024; % длина выборки - число элементов в массиве данных
% размером 2 в степени n - целое число
K = N/2; % число спектральных линий - должно быть N/2, N/4 и т.д.
% - определяет ширину спектра
Delta_t = 0.1 %sec интервал времени, в течение которого
% получаем исходный массив данных
Td = Delta_t/N; % период дискретизации входного сигнала -
% интервал через который происходит оцифровка
% сигнала в АЦП (в данном примере задана через длину выборки и длительность записи сигнала)
f_signal = 1000; %Hz частота полезного сигнала
f_p1 = 2000; %Hz частота первой помехи
f_p2 = 4000; %Hz частота второй помехи
T_signal = 1/f_signal; %sec период полезного сигнала
T_p1 = 1/f_p1; %sec период первой помехи
T_p2 = 1/f_p2; %sec период второй помехи
A_signal = 2; %V амплитуда сигнала
A_p1 = 1; %V амплитуда первой помехи
A_p2 = 0.5; %V амплитуда второй помехи
A_nois = 0.7; %V амплитуда белого шума
A_max = A_signal + A_p1 + A_p2 + A_nois; %V
% максимально возможная амплитуда результирующего сигнала
for i = 1:N, % генерация дискретного сигнала
t(i) = Delta_t*i/1024; % массив моментов выборки сигнала
% заполнение массива данных полезного сигнала
Signal(i) = A_signal*sin(2*pi/T_signal*t(i));
% заполнение массива данных первой помехи
Pomeha_1(i) = A_p1*sin(2*pi/T_p1*t(i));
% заполнение массива данных второй помехи
Pomeha_2(i) = A_p2*sin(2*pi/T_p2*t(i));
Nois(i) = A_nois*rand(1); % заполнение массива данных белого шума
% заполнение массива данных результирующего сигнала
X(i) = Signal(i) + Pomeha_1(i) + Pomeha_2(i) + Nois(i);
end
%--------------------------------------------------------------------
for i = 1:N/4, % выделение небольшой части сигнала для вывода графика
t1(i) = t(i); % массив моментов времени
X1(i) = X(i); % массив значений части сигнала
end
figure; % вывод графиков в отдельном окне
% вывод первого из двух графиков: исходный сигнал
subplot(211), plot(t1,X1,'k');
% вывод текста по оси абсцисс
xlabel('Input signal: time, sec');
ylabel('X(t)'); % вывод текста по оси ординат
% диапазон значений вывода по оси ординат
ylim([-(fix(A_max)+0.5) (fix(A_max)+0.5)]);
%--- Прямое дискретное преобразование Фурье входного сигнала ------
% цикл по числу спектральных линий (для действительной части пектра)
for i = 1:K,
S1(i) = 0; % обнуление массива действительной части спектра
S2(i) = 0; % обнуление массива мнимой части спектра
for j = 1:N, % цикл по числу элементов массива выборки
% расчет действительной части спектра
S1(i) = S1(i) + X(j)*cos(-2*pi/N*i*j);
% расчет мнимой части спектра
S2(i) = S2(i) + X(j)*sin(-2*pi/N*i*j);
end
S1(i) = S1(i)/K; % нормировка спектра
S2(i) = S2(i)/K; % нормировка спектра
% расчет модуля спектральной плотности
S(i) = sqrt(S1(i)*S1(i) + S2(i)*S2(i));
f(i) = i/N/Td; % расчет частот спектра
end
%------- Расчет критического уровня для спектральной плотности ----------
% для размера выборки N = 1024 и величины ошибки q = 0.05 достоверность
% P = 1-q = 95% порог обнаружения сигнала по критерию Стьюдента
Z1 = 28.2;
% обнуляем сумму квадратов значений дискретного входного сигнала
Sum_X = 0;
for i = 1:N, % цикл по объему выборки
% расчет суммы квадратов входного дискретного сигнала
Sum_X = Sum_X + X(i)^2;
end
Sigma_02 = Sum_X/(N-1) % расчет дисперсии входного дискретного сигнала
% расчет критического уровня, выше которого спектральная плотность значима
Critical_Level_1 = Sigma_02*Z1/N
%--------------------------------------------------------------------
%вывод второго графика: спектральной плотности и линии критического уровня
subplot(212), plot(f,S,'r',f,Critical_Level_1,'b');
% вывод текста по оси абсцисс
xlabel('Spectrum of input signal: frequency, Hz');
% вывод текста по оси ординат
ylabel('S(f)');
% ограничение диапазона значений вывода данных по оси абсцисс
xlim([0 10*K]);
% ограничение диапазона значений вывода данных по оси ординат
ylim([0 fix(A_signal)]);
%-------------- Фильтрация сигнала полосовым фильтром ---------------
f_low = 500; %Hz нижняя частота среза фильтра
f_high = 1500; %Hz верхняя частота среза фильтра
for i = 1:K, % цикл по числу спектральных линий
% если частота спектра не попадает в полосу пропускания
if (f(i)<f_low)||(f(i)>f_high),
% то обнуляем действительную и мнимую части (i)- частоты спектра
S1(i) = 0;
S2(i) = 0;
S(i) = 0;
end
end
%--------------------------------------------------------------------
figure; % вывод графиков в отдельном окне
%вывод второго графика:спектральной плотности отфильтрованного сигнала
subplot(211), plot(f,S,'r');
% вывод текста по оси абсцисс
xlabel('Spectrum of filtered signal: frequency, Hz');
% вывод текста по оси ординат
ylabel('SF(f)');
% ограничение диапазона значений вывода данных по оси абсцисс
xlim([0 10*K]);
% ограничение диапазона значений вывода данных по оси ординат
ylim([0 fix(A_signal)]);
%--- Обратное дискретное преобразование Фурье отфильтрованного сигнала
for i = 1:N, % цикл по числу выборки
Y(i) = 0; % обнуление элементов массива отфильтрованного сигнала
t(i) = i*Td; % заполнение массива моментов времени
for j = 1:K, % цикл по числу спектральных линий
Y(i) = Y(i) + S1(j)*cos(2*pi*f(j)*t(i)) + S2(j)*sin(2*pi*f(j)*t(i));
end
end
%--------------------------------------------------------------------
% выделение небольшой части отфильтрованного сигнала для вывода графика
for i = 1:N/4,
t1(i) = t(i); % массив моментов времени
Y1(i) = Y(i); % массив значений части отфильтрованного сигнала
end
% вывод второго графика в окне: отфильтрованный сигнал
subplot(212), plot(t1,Y1,'b');
% вывод текста по оси абсцисс
xlabel('Filtered signal: time, sec');
% вывод текста по оси ординат
ylabel('Y(t)');
% ограничение диапазона значений вывода данных по оси ординат
ylim([-(fix(A_max)+0.5) (fix(A_max)+0.5)]);
%--------------------------------------------------------------------
% Расчет корреляционной функции отфильтрованного сигнала с полезным сигналом
C = zeros(1,N); % обнуление массива значений корреляционной функции
for i = 1:K, % цикл по фрагменту массива отфильтрованного сигнала
for j = 1:K, % цикл по массиву полезного сигнала
C(i) = C(i) + X(i+j)*Y(j);
end
C(i) = C(i)/N; % нормировка корреляционной функции
end
% цикл восстановления корреляционной функции до длины выборки
for i = K+1:N,
C(i) = C(i-K);
end
%--------------------------------------------------------------------
figure; % вывод следующего окна для графиков
% задание интервала изменения аргумента i при выводе коррелограммы
i = 1:N;
% вывод первого графика - корреляционной функции
subplot(211), plot(i,C,'r');
% вывод текста по оси абсцисс
xlabel('Correlation function: number of count');
% вывод текста по оси ординат
ylabel('C(i)');
% ограничение диапазона значений вывода данных по оси абсцисс
xlim([1 N/4]);
% ограничение диапазона значений вывода данных по оси ординат
ylim([-1.5 1.5]);
%--- Прямое дискретное преобразование Фурье корреляционной функции ----
% цикл по числу спектральных линий (для действительной части спектра)
for i = 1:K,
S1(i) = 0;% обнуление массива действительной части спектра
S2(i) = 0;% обнуление массива мнимой части спектра
% цикл по числу элементов массива данных корреляционной функции
for j = 1:N,
% расчет действительной части спектра
S1(i) = S1(i) + C(j)*cos(-2*pi/N*i*j);
% расчет мнимой части спектра
S2(i) = S2(i) + C(j)*sin(-2*pi/N*i*j);
end
S1(i) = S1(i)/K; % нормировка спектра
S2(i) = S2(i)/K; % нормировка спектра
% расчет модуля спектральной плотности коррелограммы
S(i) = sqrt(S1(i)*S1(i) + S2(i)*S2(i));
f(i) = i/N/Td; % расчет частот коррелограммы
end
%----- Расчет критического уровня для спектральной плотности
% корреляционной функции ------------------------------------------------
% для размера выборки N = 1024 и величины ошибки q = 0.05 достоверность
% P = 1-q = 95%
Z1 = 28.2; % порог обнаружения сигнала по критерию Стьюдента
Sum_Signal = 0; % обнуление переменной суммы значений полезного сигнала
Sum_Y = 0; % обнуление переменной суммы значений отфильтрованного сигнала
for i = 1:N, % цикл по объему выборки
% расчет суммы квадратов значений полезного сигнала
Sum_Signal = Sum_Signal + Signal(i)^2;
% расчет суммы квадратов значений отфильтрованного сигнала
Sum_Y = Sum_Y + Y(i)^2;
end
Sigma_12 = Sum_Signal/(N-1) % расчет дисперсии полезного сигнала
Sigma_22 = Sum_Y/(N-1) % расчет дисперсии отфильтрованного сигнала
% расчет критического уровня выше которого
Critical_Level_2 = Sigma_12*Sigma_22*Z1/N^2
% спектральная плотность коррелограммы значима
%--------------------------------------------------------------------
%вывод второго графика в окне: спектр коррелограммы и критический уровнь
subplot(212), plot(f,S,'b',f,Critical_Level_2,'r');
% вывод текста по оси абсцисс
xlabel('Correlation spectrum: frequency, Hz');
% вывод текста по оси ординат
ylabel('SC(f)');
% ограничение диапазона значений вывода данных по оси абсцисс
xlim([0 10*K]);
% ограничение диапазона значений вывода данных по оси ординат
ylim([0 1.5]);
%--------------------------------------------------------------------
%----- расчет выборочного коэффициента корреляции исходного полезного
% сигнала с результирующим исходным сигналом-------------------------
% обнуление сумм массивов полезного и результирующего сигналов
Sum_Signal = 0;
Sum_X = 0;
for i = 1:N,
% сумма массива полезного сигнала
Sum_Signal = Sum_Signal + Signal(i);
% сумма массива отфильтрованого сигнала
Sum_X = Sum_X + X(i);
end
% расчет среднего значения полезного сигнала
m_Signal = Sum_Signal/N;
% расчет среднего значения отфильтрованного полезного сигнала
m_X = Sum_X/N;
% обнуление сумм массивов
Sum1 = 0;
Sum2 = 0;
for i = 1:N,
Sum1 = Sum1 + (Signal(i)-m_Signal)^2;
Sum2 = Sum2 + (X(i)-m_X)^2;
end
% расчет дисперсии полезного и отфильтрованного сигналов
D_Signal = Sum1/(N-1);
D_X = Sum2/(N-1);
% расчет второго смешанного центрального момента
% обнуление суммы
Sum3 = 0;
for i = 1:N,
Sum3 = Sum3 + (Signal(i)-m_Signal)*(X(i)-m_X);
end
M = Sum3/(N-1);
% расчет выборочного коеффициента корpеляции
R0 = M/sqrt(D_Signal*D_X)
%---------------------------------------------------------------
%----------- расчет выборочного коэффициента корреляции полезного и
% отфильтрованного сигналов ------------------------------------------
% обнуление сумм массивов полезного и фильтрованного сигналов
Sum_Signal = 0;
Sum_Y = 0;
for i = 1:N,
% сумма массива полезного сигнала
Sum_Signal = Sum_Signal + Signal(i);
% сумма массива отфильтрованного сигнала
Sum_Y = Sum_Y + Y(i);
end
% расчет среднего значения полезного сигнала
m_Signal = Sum_Signal/N
% расчет среднего значения отфильтрованного полезного сигнала
m_Y = Sum_Y/N
% обнуление сумм массивов
Sum1 = 0;
Sum2 = 0;
for i = 1:N,
Sum1 = Sum1 + (Signal(i)-m_Signal)^2;
Sum2 = Sum2 + (Y(i)-m_Y)^2;
end
% расчет дисперсии полезного и отфильтрованного сигналов
D_Signal = Sum1/(N-1)
D_Y = Sum2/(N-1)
% расчет второго смешанного центрального момента
% обнуление суммы
Sum3 = 0;
for i = 1:N,
Sum3 = Sum3 + (Signal(i)-m_Signal)*(Y(i)-m_Y);
end
M = Sum3/(N-1)
% расчет выборочного коэффициента корреляции
R = M/sqrt(D_Signal*D_Y)
%---------------------------------------------------------------
Фрагмент программы моделирования цифровой фильтрации зашумленного амплитудно-модулированного сигнала:
% Фильтрация амплитудно-модулированного сигнала
clc;% очистить командное окно
clear all;% очистить рабочее пространство
digits(4);% точность расчетов 4 знака после запятой
%--------------------------------------------------------------------
N = 1024;% длина выборки - число элементов в массиве данных
% размером 2 в степени n - целое число
K = N/2;% число спектральных линий - должно быть N/2, N/4 и т.д.
% - определяет ширину спектра
Delta_t = 0.1 %sec интервал времени, в течение которого
% получаем исходный массив данных
Td = Delta_t/N; % период дискретизации входного сигнала -
% интервал через который происходит оцифровка
% сигнала в АЦП (в данном примере задана через длину выборки и длительность
% записи сигнала)
f_mod = 50; %Hz частота амплитудной модуляции
f_signal = 1000; %Hz частота полезного сигнала
f_p1 = 2000; %Hz частота первой помехи
f_p2 = 4000; %Hz частота второй помехи
T_mod = 1/f_mod; %sec период амплитудной модуляции
T_signal = 1/f_signal; %sec период полезного сигнала
T_p1 = 1/f_p1; %sec период первой помехи
T_p2 = 1/f_p2; %sec период второй помехи
A_signal = 2; %V максимальная амплитуда полезного сигнала
A_p1 = 1; %V амплитуда первой помехи
A_p2 = 0.5; %V амплитуда второй помехи
A_nois = 0.5; %V амплитуда белого шума
% максимально возможная амплитуда результирующего сигнала
A_max = A_signal + A_p1 + A_p2 + A_nois; %V
for i = 1:N, % генерация дискретного сигнала
t(i) = Delta_t*i/1024; % заполнение массива моментов выборки сигнала
% заполнение массива данных полезного сигнала
Signal(i) = A_signal*cos(2*pi/T_mod*t(i))*sin(2*pi/T_signal*t(i));
% заполнение массива данных первой помехи
Pomeha_1(i) = A_p1*sin(2*pi/T_p1*t(i));
% заполнение массива данных второй помехи
Pomeha_2(i) = A_p2*sin(2*pi/T_p2*t(i));
Nois(i) = A_nois*rand(1);% заполнение массива данных белого шума
% заполнение массива данных результирующего сигнала
X(i) = Signal(i) + Pomeha_1(i) + Pomeha_2(i) + Nois(i);
end
%--------------------------------------------------------------------
Фрагмент программы моделирования цифровой фильтрации зашумленного частотно-модулированного сигнала:
% Фильтрация частотно-модулированного сигнала
clc;% очистить командное окно
clear all;% очистить рабочее пространство
digits(4);% точность расчетов 4 знака после запятой
%--------------------------------------------------------------------
N = 1024; % длина выборки - число элементов в массиве данных
% размером 2 в степени n - целое число
K = N/2; % число спектральных линий - должно быть N/2, N/4 и т.д.
% - определяет ширину спектра
Delta_t = 0.1 %sec интервал времени, в течение которого
% получаем исходный массив данных
Td = Delta_t/N;% период дискретизации входного сигнала -
% интервал через который происходит оцифровка
% сигнала в АЦП (в данном примере задана через длину выборки и длительность
% записи сигнала)
f_mod = 100; %Hz частота модулирующего сигнала
f_signal = 1000; %Hz несущая частота сигнала (немодулированного сигнала)
f_p1 = 2000; %Hz частота первой помехи
f_p2 = 4000; %Hz частота второй помехи
T_mod = 1/f_mod; %sec период модулирующего сигнала
% амплитуда изменения периода колебаний несущего сигнала
Delta_T_mod = 0.0003*T_mod; %sec
T_signal_0 = 1/f_signal; %sec период немодулированного полезного сигнала
T_p1 = 1/f_p1; %sec период первой помехи
T_p2 = 1/f_p2; %sec период второй помехи
A_signal = 2; %V амплитуда полезного сигнала
A_p1 = 0.5; %V амплитуда первой помехи
A_p2 = 0.25; %V амплитуда второй помехи
A_nois = 0.25; %V амплитуда белого шума
% максимально возможная амплитуда результирующего сигнала
A_max = A_signal + A_p1 + A_p2 + A_nois; %V
for i = 1:N, % генерация дискретного сигнала
t(i) = Delta_t*i/1024; % заполнение массива моментов выборки сигнала
% заполнение массива периода модулированного полезного сигнала
T_signal(i) = T_signal_0 + Delta_T_mod*sin(2*pi/T_mod*t(i));
% заполнение массива данных полезного сигнала
Signal(i) = A_signal*sin(2*pi/T_signal(i)*t(i));
% заполнение массива данных первой помехи
Pomeha_1(i) = A_p1*sin(2*pi/T_p1*t(i));
% заполнение массива данных второй помехи
Pomeha_2(i) = A_p2*sin(2*pi/T_p2*t(i));
% заполнение массива данных белого шума
Nois(i) = A_nois*rand(1);
% заполнение массива данных результирующего сигнала
X(i) = Signal(i) + Pomeha_1(i) + Pomeha_2(i) + Nois(i);
end
%--------------------------------------------------------------------