Всем доброго времени суток и с наступившими праздниками!
В статье пойдёт речь о построении графиков в MATLAB для моей курсовой по дисциплине радиотехнические цепи и сигналы. Никаких мегакрутых и инновационных алгоритмов и тайных знаний здесь нет. Статья предназначена для помощи студентам с потока в оформлении курсовой работы. Итак, поехали:
Для тех, кто не хочет читать, сразу выложу свой код, а потом постараюсь подробно всё описать.
function PTC
%% 1 пункт
figure
x=[];
y=[];
for t=-10:0.01:10
x=[x,t];
t=t/1000;
y=[y,sig(t)];
end
plot(x,y), grid
xlabel('t, мс');
ylabel('s(t), В');
title('Сигнал');
%% Спектры
figure
x=[];
y=[];
subplot(2,1,1)
for w=-15:0.01:15
x=[x,w];
w=w*1000;
y=[y,1000*S(w)];
end
plot(x,y), grid
xlabel('\omega, рад/мс');
ylabel('|S(\omega)|, В/кГц');
title('Амплитудный спектр');
%%
y=[];
subplot(2,1,2)
for w=-15:0.01:15
w=w*1000;
y=[y,fi(w)];
end
plot(x,y), grid
xlabel('\omega, рад/мс');
ylabel('\phi(\omega), \circ');
title('Фазовый спектр');
%% 2 пункт
% Sп(t)
figure
x=[];
y=[];
T=30*10^-3;
for t=-80:0.1:80
s=0;
x=[x,t];
t=t/1000;
for n=-10:10
s=s+sig(t-n*T);
end
y=[y,s];
end
plot(x,y,'LineWidth' ,2),grid
xlabel('t, мс');
ylabel('s_п(t), В');
title('Временная диаграмма периодического сигнала');
%% Спектры
figure
x=[];
y=[];
subplot(2,1,1,'xlim',[0,10])
for n=0:50
w=2*pi*n/T; %ось х мс
x=[x,w/1000];
s=2/T*S(w);
if n==0
s=s/2;
end
y=[y,s];
line([w/1000 w/1000],[0 s]);
line(w/1000, s,'marker','.','markersize' ,10)
end
hold on
plot(x,y,'--','Color',[1 0 0]);
grid,
xlabel('\omega, рад/мс');
ylabel('\{A_n\}, В');
title('Амплитудный спектр');
%%
subplot(2,1,2,'xlim',[0,10])
for n=0:50
w=2*pi*n/T;
line([w/1000 w/1000],[0 fi(w)]);
line(w/1000, fi(w),'marker','.','markersize' ,10);
end
x=[];
y=[];
for w=0:0.01:10
x=[x,w];
w=w*1000;
y=[y,fi(w)];
end
hold on
plot(x,y,'--','Color',[1 0 0]);
grid
xlabel('\omega, рад/мс');
ylabel('\{\phi_n\}, \circ');
title('Фазовый спектр');
%% Апроксимация n=10
figure
x=[];
y=[];
subplot(3,1,1)
for t=-15:0.01:15
x=[x,t];
t=t/1000;
s=2/15;
for n=1:10
w=2*pi*n/T; %ось х мс
s=s+2/T*S(w)*cos(w*t+pi/180*fi(w));
end
y=[y,s];
end
plot(x,y)
grid
xlabel('t, мс');
ylabel('s_п(t), В');
title('n=10');
%% n=25
x=[];
y=[];
subplot(3,1,2)
for t=-15:0.01:15
x=[x,t];
t=t/1000;
s=2/15;
for n=1:25
w=2*pi*n/T; %ось х мс
s=s+2/T*S(w)*cos(w*t+pi/180*fi(w));
end
y=[y,s];
end
plot(x,y)
grid
xlabel('t, мс');
ylabel('s_п(t), В');
title('n=25');
%% n=50
x=[];
y=[];
subplot(3,1,3)
for t=-15:0.01:15
x=[x,t];
t=t/1000;
s=2/15;
for n=1:50
w=2*pi*n/T; %ось х мс
s=s+2/T*S(w)*cos(w*t+pi/180*fi(w));
end
y=[y,s];
end
plot(x,y)
grid
xlabel('t, мс');
ylabel('s_п(t), В');
title('n=50');
%% 3 пункт
figure
x=[];
y=[];
v=[];
w=20000;
for t=-5:0.001:5
x=[x,t];
t=t/1000;
y=[y,sig(t)*cos(w*t)];
v=[v,sig(t)];
end
hold on;
plot(x,y),grid
plot(x,v,'--','Color',[1 0 0])
plot(x,-v,'--','Color',[1 0 0])
xlabel('t, мс');
ylabel('u(t), В');
title('Временная диаграмма непериодического радиосигнала');
%% положительные частоты
figure
x=[];
y=[];
subplot(2,1,1)
for w=490:0.001:510
x=[x,w];
w=w*1000;
y=[y,500*S(w-500000)];
end
plot(x,y), grid
xlabel('\omega, рад/мс');
ylabel('|U(\omega)|, В/кГц');
title('Амплитудный спектр');
%%
y=[];
subplot(2,1,2)
for w=490:0.001:510
w=w*1000;
y=[y,fi(w-500000)];
end
plot(x,y), grid
xlabel('\omega, рад/мс');
ylabel('\phi_U(\omega), \circ');
title('Фазовый спектр');
%% отрицательные частоты
figure
x=[];
y=[];
subplot(2,1,1)
for w=-510:0.001:-490
x=[x,w];
w=w*1000;
y=[y,500*S(w+500000)];
end
plot(x,y), grid
xlabel('\omega, рад/мс');
ylabel('|U(\omega)|, В/кГц');
title('Амплитудный спектр');
%%
y=[];
subplot(2,1,2)
for w=-510:0.001:-490
w=w*1000;
y=[y,fi(w+500000)];
end
plot(x,y), grid
xlabel('\omega, рад/мс');
ylabel('\phi_U(\omega), \circ');
title('Фазовый спектр');
%% 4 пункт
% Sп(t)
figure
x=[];
y=[];
v=[];
s1=0;
s2=0;
w0=15000;
for t=-40:0.001:40
x=[x,t];
t=t/1000;
for n=-1:1
s1=s1+sig(t-n*T)*cos(w0*t);
s2=s2+sig(t-n*T);
end
y=[y,s1];
v=[v,s2];
s1=0;
s2=0;
end
hold on
plot(x,y),grid
plot(x,v,'--','Color',[1 0 0])
plot(x,-v,'--','Color',[1 0 0])
xlabel('t, мс');
ylabel('u_п(t), В');
title('Временная диаграмма периодического радиосигнала');
%% Спектры пачки радиоимпульсов
figure
x=[];
y=[];
subplot(2,1,1,'xlim',[490 510])
for n=-50:50
w=2*pi*n/T; %ось х мс
x=[x,w/1000+500];
y=[y,1/T*S(w)];
line([w/1000+500 w/1000+500],[0 1/T*S(w)]);
line(w/1000+500, 1/T*S(w),'marker','.','markersize' ,10)
end
hold on
plot(x,y,'--','Color',[1 0 0]);
grid,
xlabel('\omega, рад/мс');
ylabel('\{V_n\}, В');
title('Амплитудный спектр');
%%
x=[];
y=[];
subplot(2,1,2,'xlim',[490 510])
for n=-50:50
w=2*pi*n/T; %ось х мс
x=[x,w/1000+500];
y=[y,fi(w)];
line([w/1000+500 w/1000+500],[0 fi(w)]);
line(w/1000+500, fi(w),'marker','.','markersize' ,10);
end
x=[];
y=[];
for w=-10:0.01:10
x=[x,w+500];
w=w*1000;
y=[y,fi(w)];
end
hold on
plot(x,y,'--','Color',[1 0 0]);
grid
xlabel('\omega, рад/мс');
ylabel('\{\psi_n\}, \circ');
title('Фазовый спектр');
%% 5 пункт
%АКФ
figure
x=[];
y=[];
for t=-10:0.01:10
x=[x,t];
t=t/1000;
y=[y,1000*R(t)];
end
plot(x,y), grid
xlabel('\tau, мс');
ylabel('R_s(\tau), В^2мc');
title('Автокорелляционная функция сигнала');
%% АКФ радиоимпульса
figure
x=[];
y=[];
for t=-10:0.01:10
x=[x,t];
t=t/1000;
y=[y,1/2*1000*R(t)*cos(w0*t)];
end
hold on;
plot(x,y), grid
x=[];
y=[];
for t=-10:0.01:10
x=[x,t];
t=t/1000;
y=[y,1/2*1000*R(t)];
end
plot(x,y,'--','Color',[1 0 0]);
plot(x,-y,'--','Color',[1 0 0]);
xlabel('\tau, мс');
ylabel('R_u(\tau), В^2мc');
title('Автокорелляционная функция радиосигнала');
%% 6 пункт
% ачх фчх
figure
x=[];
y=[];
subplot(2,1,1)
for wt=-5:0.01:5
x=[x,wt];
y=[y,(1/36*wt^2/(1+wt^2))^0.5];
end
plot(x,y), grid
xlabel('\omega\tau, рад');
ylabel('|H(\omega)|');
title('АЧХ');
%%
y=[];
subplot(2,1,2)
for wt=-5:0.01:5
y=[y,(180/pi)*(pi/2*sign(wt)-atan(wt))];
end
plot(x,y), grid
xlabel('\omega\tau, рад');
ylabel('\phi_H(\omega), \circ');
title('ФЧХ');
%% ИХ ПХ g1
figure
x=[];
y=[];
subplot(3,1,1)
for t=0:0.001:1
x=[x,t];
t=t/1000;
y=[y,(-1/(6*0.0002)*sigma(t)*exp(-t/0.0002))/1000];
end
plot(x,y,'LineWidth' ,2)
grid
xlabel('t, мс');
ylabel('h(t), мс^-1');
title('Импульсная характеристика');
%% ПХ
x=[];
y=[];
subplot(3,1,2)
for t=0:0.001:1
x=[x,t];
t=t/1000;
y=[y,g(t)];
end
plot(x,y,'LineWidth' ,2)
grid
xlabel('t, мс');
ylabel('g(t)');
title('Переходная характеристика');
%% Реакция на лин нар возд
x=[];
y=[];
subplot(3,1,3)
for t=0:0.001:1
x=[x,t];
t=t/1000;
y=[y,(1/6*0.0002*(1-exp(-t/0.0002)))];
end
plot(x,y,'LineWidth' ,2)
grid
xlabel('t, мс');
ylabel('g_1(t)');
title('g_1(t)');
%% Запаздывание
figure
x=[];
y=[];
for t=-2:0.001:8
x=[x,t];
t=t/1000;
y=[y,sig(t-0.002)];
end
plot(x,y,'LineWidth' ,2)
grid
xlabel('t, мс');
ylabel('s_з(t), В');
title('Запаздывающий сигнал');
%% uвх
figure
x=[];
y=[];
for t=-2:0.001:8
x=[x,t];
t=t/1000;
s=sigma(t)-sigma(t-0.002)+sigma(t-0.004)-sigma(t-0.006);
y=[y,s];
end
plot(x,y,'LineWidth' ,1)
grid
xlabel('t, мс');
ylabel('u_{вых}(t), В');
title('Cигнал на входе цепи');
%% uвых
figure
x=[];
y=[];
for t=0:0.001:8
x=[x,t];
t=t/1000;
s=g(t)-g(t-0.002)+g(t-0.004)-g(t-0.006);
y=[y,s];
end
plot(x,y,'LineWidth' ,1)
grid
xlabel('t, мс');
ylabel('u_{вых}(t), В');
title('Cигнал на выходе цепи');
end
%%
function f=S(w)
t1=0.002;
f=2*t1*abs(sinc(w*t1/2)*cos(w*t1));
end
%%
function f=fi(w)
t1=0.002;
f=(180/pi)*(pi/2*(1-sign(sinc(w*t1/2)*cos(w*t1)))*sign(w)-w*t1/2);
end
%%
function f=sig(t)
t1=0.002;
f=(rect((t+t1/2)/t1)+rect((t-3*t1/2)/t1));
end
%%
function f=R(t)
t1=0.002;
f=(2*t1*(1-abs(t)/t1)*rect(t/(2*t1))+t1*(1-abs(abs(t)-2*t1)/t1)*rect((abs(t)-2*t1)/(2*t1)));
end
%%
function f=g(t)
f=1/6*sigma(t)-1/6*sigma(t)*(1-exp(-t/0.0002));
end
%%
function S=sigma(x)
if x<0
S=0;
else
S=1;
end
end
%%
function S=rect(x)
if abs(x)>0.5
S=0;
end
if abs(x)<0.5
S=1;
end
if abs(x)==0.5
S=0.5;
end
end
%%
function S=sinc(x)
if x==0
S=1;
else
S=sin(x)/x;
end
end
%%
function S=sign(x)
if x<0
S=-1;
end
if x>0
S=1;
end
if x==0
S=0;
end
end
Введение
Итак, вы установили MATLAB. Теперь его надо запустить.
После этого в поле Command Window нужно написать edit, чтобы открыть редактор, где и будем писать программу.
Расположение окон может быть несколько иным, но суть не меняется. Теперь будем разбираться с самой программой.
%% 1 пункт
figure
x=[];
y=[];
for t=-10:0.01:10
x=[x,t];
t=t/1000;
y=[y,sig(t)];
end
plot(x,y), grid
xlabel('t, мс');
ylabel('s(t), В');
title('Сигнал');
figure нужен для создания окна, в котором и будут построены нужные графики.
Далее создаём пустые массивы x и y, которые будут содержать координаты точек.
Теперь нужно заполнить эти массивы. Для этого создадим цикл for. Конструкция t=-10:0.01:10 означает: переменная t будет меняться от -10 до 10 с шагом 0.01(по умолчанию шаг равен 1). Чем меньше шаг, тем более плавным получается график, но медленнее строится, т.к. операций становится больше.
x=[x,t] – производит добавление элемента к массиву с сохранением предыдущих. В данном примере это выглядит так:
Изначально массив пустой x=[];
1. Первая итерация - x=[-10] ;
2. x=[-10 -9,99];
3. x=[-10 -9,99 -9,98];
…
n. x=[-10 -9,99 -9,98 … 9,99 10];
t=t/1000 – рисовать будем не в секундах, а в миллисекундах. Можно этого не делать, тогда в условие цикла будет такое: t=-0.010:0.00001:0.010. От нулей рябит в глазах, плюс выше вероятность ошибки.
y=[y,sig(t)]; - заполняем массив y. sig(t) – функция, которую мы определим в конце. Она будет возвращать значение функции заданного сигнала, т.е t – это x, а sig(t) – y.
Вот как я описал sig(t):
function f=sig(t) t1=0.002; f=(rect((t+t1/2)/t1)+rect((t-3*t1/2)/t1)); end
Аналогично ввожу функцию rect
function S=rect(x)
if abs(x)>0.5
S=0;
end
if abs(x)<0.5
S=1;
end
if abs(x)==0.5
S=0.5;
end
end
abs(x) возвращает модуль числа x.
Команда plot(x, y) соответствует построению функции, когда одномерный массив x соответствует значениям аргумента, а одномерный массив y - значениям функции.
grid добавляет на график сетку
xlabel и уlabel – добавление подписей осей.
title – название графика
Результат работы этой части программы(нажать F5 для запуска)
Можно сохранять картинку. Жмём на соответствующую иконку, выбираем формат bmp или jpg, сохраняем.
Аналогично рисуем спектры.
figure
x=[];
y=[];
subplot(2,1,1)
for w=-15:0.01:15
x=[x,w];
w=w*1000;
y=[y,1000*S(w)];
end
plot(x,y), grid
xlabel('\omega, рад/мс');
ylabel('|S(\omega)|, В/кГц');
title('Амплитудный спектр');
%%
y=[];
subplot(2,1,2)
for w=-15:0.01:15
w=w*1000;
y=[y,fi(w)];
end
plot(x,y), grid
xlabel('\omega, рад/мс');
ylabel('\phi(\omega), \circ');
title('Фазовый спектр');
subplot(2,1,1) и subplot(2,1,2) позволяют в одном окне рисовать графики в первой и второй половине окна.
w=w*1000; - как с секундами, только в одном рад/с тысяча рад/мс
По тем же причинам перемножаем значения амплитудного спектра на 1000.
\omega – так можно напечатать символ омега малое в матлабе.
Нужно отметить, чтобы нарисовать фазовый спектр в градусах, нужно домножить всё выражение на 180/3.1415.
Как это выглядит
Переходим к следующему пункту. Далее я не буду приводить целые куски кода, поэтому советую перенести код из начала статьи в матлаб.
T – период повторения импульса.
for n=-10:10
s=s+sig(t-n*T);
end
Этот цикл позволяет получить выражение для пачки импульсов из выражения непериодического сигнала. В курсовой n изменяется от минус до плюс бесконечности. Но для построения нам будет предостаточно 20 периодов(лишние всё равно не будут видны, границы графика определяется циклов for выше. Я подгонял эти границы, чтобы увидеть 5 периодов).
При построении спектральных диаграмм пришлось по-другому ограничивать эту область ( subplot(2,1,1,'xlim',[0,10]) )
Спектр периодического сигнала дискретен, для получения выражений нужно использовать связь непериодического и периодического сигналов. Для построения спектральных линий буду использовать оператор line. n в цикле for – количество этих линий. Подбирается вручную.
line([x1 x2],[y1 y2])
line([w/1000 w/1000],[0 s]) – рисует перпендикуляр от оси абсцисс.
line(w/1000, s,'marker','.','markersize' ,10) - изображает точку('marker','.'), 'markersize' ,10 –толщина точки.
hold on – позволяет строить несколько функций на одном графике.
plot(x,y,'--','Color',[1 0 0]); - рисует пунктирную линию («--»), красным цветом('Color',[r g b]).
При построении графиков нужно убедиться, что фаза выражена в радианах.
В третьем и четвёртом пунктах проводится исследования радиосигналов
w0=15000; w=20000 – частота несущей, которая будет использоваться при построении(подбирается, чтобы было наглядно)
500000 – реальная частота несущей(находим в курсовой)
Далее сложностей уже не должно возникнуть.
Удачи при выполнении курсовой!
Леонид Делов, 3 января 2019
lad_kursk@mail.ru
https://vk.com/l_deloff