Как построить ачх в матлабе

Matlab (Теория автоматического управления-Control System Toolbox)

Розрахунок одноконтурної системи автоматичного керування за допомогою пакету Matlab

Текст программы :

clc
clear
sys=tf(1.2,[85,1]);
sys.ioDelay=35;
get(sys);

figure(1)

subplot(2,1,1)
step(sys),grid on

subplot(2,1,2)
impulse(sys),grid on

figure(2)
[Re_,Im_,Omega]=nyquist(sys,[0:0.001:35]);
Re=[];
Im=[];
Re(:,1,1)=Re_(1,: );
Im(:,1,1)=Im_(1,: );
plot(Re,Im,Re,Im),grid on

figure(3)
[Amp,Phasep,Omega]=bode(sys,[0:0.001:0.1]);
Am=[];
Phase=[];
Am(:,1)=Amp(1,1,: );
Phase(:,1)=Phasep(1,1,: );
subplot(2,1,1)
plot(Omega,Am),grid on
subplot(2,1,2)
plot(Omega,Phase),grid on

c0=0;
figure(4)
while c0>=0
A=1-(T*m*w);
B=-(T*w);
c1=((A*m-B)*sin(w*tau)-(A+m*B)*cos(w*tau))/(k*exp(m*w*tau));
c0=((A*sin(w*tau)-B*cos(w*tau))/(k*exp(m*w*tau)))*((1+m^2)*w);
vc0=[vc0,c0];
vc1=[vc1,c1];
w=w+hw;
end

n=length(vc0);
n=length(vc1);
vc0(n)=[];
vc1(n)=[];
[c0max,n]=max(vc0);
c0opt=vc0(n+1)
c1opt=vc1(n+1)
plot(vc1,vc0,vc1,vc0,’r.’,c1opt,c0opt,’green*’),grid on

figure(5)
w0=tf([k],[T,1]);
w0.ioDelay=tau;
w0ss=ss(w0);
wp=tf([c1opt,c0opt],[1 0]);
wpss=ss(wp);
w=feedback(w0ss,wpss);
t=[0:1:600];
[ht,t]=step(w,t);
plot(t,ht),grid on ,hold on

n=length(ht);
[Ddin,ndin]=max(ht)
tdin=[t(ndin),t(ndin)];
Ddinh=[0,Ddin];
plot(tdin,Ddinh,’k’,tdin(1),Ddinh(1),’kv’,tdin(2),Ddinh(2),’k^’)

s4et4ik=0;
for i=[2:n-1];
if (ht(i)>ht(i-1))&(ht(i)>ht(i+1))
s4et4ik=s4et4ik+1;
if s4et4ik==2
Dstat=ht(i)
Treg=t(i)
end
end
end
plot(Treg,Dstat,’red*’), hold on
Dstath=[-Dstat];
IntArea=sum(abs(ht)*(t(2)-t(1)))

figure (6)
logvec=logspace(-2,-0.876,2000);
nyquist(w0*wp,logvec),grid on ;

figure (7)
a0 = 0.032;
a1 = 1.698;
a2 =85;
Re=[];Im=[];
for w=0.0001:0.001:0.45,
Njw= a2*((w*j)^2)+a1*(w*j)+(a0);
Re = real(Njw);
Im = imag(Njw);
plot(Re, Im, ‘k.’)
xlabel(‘Re(W)’)
ylabel(‘Im(W)’)
hold on
end
hold off
grid on
axis([-0.45 0.45 0.45 0.45])

A=[a1,0;a2,a0];
D=det(A);
if D>0&a1>0
disp (‘система устойчива’);
else
disp (‘ система не устойчива’);
end

Поделиться22012-05-05 07:59:42

• Автор: Admin
• Администратор

Переходной характеристикой называется такая характеристика объекта (или системы), которая показывает изменение выходной величины во времени звена, объекта регулятора или системы при подаче на вход единичного скачка . Единичный скачек отвечает 100%-ому мгновенному изменению управляющего сигнала или возмущающего воздействия.

Формируем модель передаточной функции:

Строим график переходной характеристики
figure(1)
subplot(2,1,1)
step(sys),grid on

Поделиться32012-05-05 08:07:39

• Автор: Admin
• Администратор

Импульсной характеристикой или весовой функцией называют такую функцию, которая описывает выходную величину, если на вход подана функция Дирака – единичная импульсная функция – мгновенная, бесконечная по амплитуде, величина с единичной площадью.

subplot(2,1,2)
impulse(sys),grid on

Поделиться42012-05-05 08:13:39

• Автор: Admin
• Администратор

Комплексной передаточной функцией называется отношение выходной величины, преобразованной по Фурье, к входной величине, преобразованной по Фурье в режиме незатухающих гармонических колебаний.
Амплитудно-частотной характеристикой (АЧХ) называется отношение амплитуды выходного гармонического сигнала к амплитуде входного гармонического сигнала.
Фазо-частотной характеристикой (ФЧХ) называется разность между входными и выходными колебаниями.
Амплитудо-фазо-частотная характеристика (АФЧХ) представляет собой геометрическое место точек на комплексной плоскости концов вектора при изменении от 0 до бесконечности.

Построение АЧХ, ФЧХ:
figure(3)
[Amp,Phasep,Omega]=bode(sys,[0:0.001:0.1]);
Am=[];
Phase=[];
Am(:,1)=Amp(1,1,: );
Phase(:,1)=Phasep(1,1,: );
subplot(2,1,1)
plot(Omega,Am),grid on
subplot(2,1,2)
plot(Omega,Phase),grid on

Поделиться52012-05-05 08:17:54

• Автор: Admin
• Администратор

Построение АФЧХ (Matlab)

figure(2)
[Re_,Im_,Omega]=nyquist(sys,[0:0.001:35]);
Re=[];
Im=[];
Re(:,1,1)=Re_(1,: );
Im(:,1,1)=Im_(1,: );
plot(Re,Im,Re,Im),grid on

Поделиться62012-05-05 08:24:00

• Автор: Admin
• Администратор

Расчет оптимальных настроек ПИ регулятора для одноконтурной системы методом расширенных АФЧХ.
Передаточная функция ПИ регулятора описывается выражением:
где C1 и C0 соответствуют пропорциональному и интегральному коэффициентам регулятора.
В общем случае для объекта 1-ого порядка C1 и C0 имеют вид:
Далее, задача сводится к построению графика линии равной степени затухания, т.е. зависимости С0(С1) , где С0=С0(W), C1=C1(W) . График строится до тех пор, пока первый раз не пересечет ось абсцисс, т.е. пока C0i(w)>0 . Определяем координаты точки, соответствующие максимальному значению C0(w) . Следующая за ней точка с координатами по ω и будет соответствовать оптимальным параметрам ПИ регулятора.

Текст программы
m=0.366; -степень колебательности (выбирается в зависимости от степени затухания)
tau=35; -запаздывание
k=1.2;
T=85; — постоянная времени

figure(5)
while c0>=0
A=1-(T*m*w);
B=-(T*w);
c1=((A*m-B)*sin(w*tau)-(A+m*B)*cos(w*tau))/(k*exp(m*w*tau));
c0=((A*sin(w*tau)-B*cos(w*tau))/(k*exp(m*w*tau)))*((1+m^2)*w);
vc0=[vc0,c0];
vc1=[vc1,c1];
w=w+hw;
end

n=length(vc0);
n=length(vc1);
vc0(n)=[];
vc1(n)=[];
[c0max,n]=max(vc0);
c0opt=vc0(n+1)
c1opt=vc1(n+1)
plot(vc1,vc0,vc1,vc0,’r.’,c1opt,c0opt,’green*’),grid on

График:

Поделиться72012-05-05 08:41:54

• Автор: Admin
• Администратор

Критерий Рауса-Гурвица
Из теории известно, что критерий сводится к составлению определителя Гурвица и определению положительности соответствующих миноров. Для АСР 1-ого и 2-ого порядка с ПИ регулятором (как в прямом контуре, так и в контуре обратной связи) характеристическое уравнение будет иметь вид

Если поставить коэф-ты в соответствии с записью векторов в MATLAB, то к соответствующему индексу надо прибавить 1.
М.б. два варианта решения проблемы:
1. Из коэф-тов вектора составляется определитель в виде квадратной матрицы и вычисляются определители всех диагональных алгебраических дополнений, путем пошагового удаления последней строки и последнего столбца. Все они должны быть больше 0 (положительны).
Длина вектора определяется функцией length(Nvec), а определитель вычисляется функцией det(mat).

Прямое сравнение коэффициентов. Для системы 1-ого порядка:a0>0,a1>0,a2>0.

Для системы 2-ого порядка: a0>0,a1>0,(a1*a2)-(a0*a3)>0, a3>0.
Откуда можно сделать вывод, что система управления устойчива/

A=[a1,0;a2,a0];
D=det(A);
if D>0 & a1>0
disp (‘система устойчива’);
else
disp (‘ система не устойчива’);
end

Поделиться82012-05-05 08:47:33

• Автор: Admin
• Администратор

Критерий Михайлова
Для того чтобы система была устойчива, необходимо и достаточно, чтобы при ее кривая Михайлова, начинаясь с положительной вещественной полуплоскости, последовательно обходила n квадрантов в положительном направлении (против часовой стрелки). Следует отметить, что кривые Михайлова устойчивых систем не пересекают начало координат и уходят на бесконечность в n-ом квадранте.
Для использования этого критерия в характеристическом уравнении оператор Лапласа заменяется на jw , Вследствие чего получим уравнение годографа Михайлова. Замена оператора Лапласа в системе MATLAB осуществляется с помощью цикла for

Выделяя вещественную и мнимую части функциями real и imag, строим годограф Михайлова в координатах .

Код:
figure (7)
a0 = 0.032;
a1 = 1.698;
a2 =85;
Re=[];Im=[];
for w=0.0001:0.001:0.45,
Njw= a2*((w*j)^2)+a1*(w*j)+(a0);
Re = real(Njw);
Im = imag(Njw);
plot(Re, Im, ‘k.’)
xlabel(‘Re(W)’)
ylabel(‘Im(W)’)
hold on
end
hold off
grid on
axis([-0.45 0.45 0.45 0.45])

Поделиться92012-05-05 08:50:50

• Автор: Admin
• Администратор

Критерий Найквиста
Критерий можно записать так: если кривая АФЧХ разомкнутой АСР не охватывает точку (–1; j0), то замкнутая АСР является устойчивой.
Поэтому, строим диаграмму Найквиста для системы с помощью функции Nyquist
Код:
figure (6)
logvec=logspace(-2,-0.876,2000);
nyquist(w0*wp,logvec),grid on ;

Поделиться102012-05-05 09:01:21

• Автор: Admin
• Администратор

Построение переходного процесса в АСР (Способ не работает в версиях Matlab 6.5 и ниже)
Качество переходного процесса
Переходной процесс – это изменение регулируемой величины во времени при переходе из одного установившегося состояния в другое.
Характер переходного процесса системы зависит от вида возмущающего воздействия и начальных условий. Для того, чтобы можно было сравнить между собой системы по характеру переходного процесса, из возможных воздействий выбирают типовые или наиболее неблагоприятные и определяют кривую переходного процесса при нулевых начальных условиях. В качестве типовых воздействий при анализе динамики АСР обычно используют единичное ступенчатое воздействие, единичный импульс, линейно возрастающее воздействие, синусоидальное воздействие и др. Для большинства систем типовым и наиболее неблагоприятным является воздействие вида единичной ступенчатой функции . Реакция системы на единичное ступенчатое воздействие при нулевых начальных условиях называется переходной функцией системы. Переходная функция системы оценивается с помощью совокупности характеристик, называемых показателями качества переходного процесса. Обычно различают следующие показатели качества переходного процесса системы:
1. Динамическая ошибка регулирования – максимальное отклонение регулируемой величины. Не должна превышать наперед заданного значения динамической ошибки.
2. Статическая ошибка – это остаточное отклонение регулируемой величины от заданного значения. Не должна превышать наперед заданного значения.
3. Время регулирования или время переходного процесса – это время от начального изменения выходной величины до ее нового установившегося значения. Время регулирования характеризует быстроту затухания переходного процесса. В этом случае требуется минимизировать.

figure(5)
w0=tf([k],[T,1]);
w0.ioDelay=tau;
w0ss=ss(w0);
wp=tf([c1opt,c0opt],[1 0]);
wpss=ss(wp);
w=feedback(w0ss,wpss);
t=[0:1:600];
[ht,t]=step(w,t);
plot(t,ht),grid on,hold on

n=length(ht);
[Ddin,ndin]=max(ht)
tdin=[t(ndin),t(ndin)];
Ddinh=[0,Ddin];
plot(tdin,Ddinh,’k’,tdin(1),Ddinh(1),’kv’,tdin(2),Ddinh(2),’k^’)

s4et4ik=0;
for i=[2:n-1];
if (ht(i)>ht(i-1))&(ht(i)>ht(i+1))
s4et4ik=s4et4ik+1;
if s4et4ik==2
Dstat=ht(i)
Treg=t(i)
end
end
end
plot(Treg,Dstat,’red*’), hold on
Dstath=[-Dstat];
IntArea=sum(abs(ht)*(t(2)-t(1))) (подинтегральная площадь графика)

3. Частотные характеристики систем автоматического управления (АФЧХ, ЛАХ, ФЧХ) ч. 3.1

Лекции по курсу «Управление Техническими Системами» читает Козлов Олег Степанович на кафедре «Ядерные реакторы и энергетические установки» факультета «Энергомашиностроения» МГТУ им. Н.Э. Баумана. За что ему огромная благодарность!

Данные лекции готовятся к публикации в виде книги, а поскольку здесь есть специалисты по ТАУ, студенты и просто интересующиеся предметом, то любая критика приветствуется.

В этом разделе мы будем изучать частотные характеристики. Тема сегодняшней статьи:
3.1. Амплитудно-фазовая частотная характеристика: годограф, АФЧХ, ЛАХ, ФЧХ

Будет интересно, познавательно и жестко.

3.1. Амплитудно-фазовая частотная характеристика: годограф АФЧХ, ЛАХ, ФЧХ

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

Формула синусоидального воздействия может быть записана как:

— сдвиг фазы (нередко называют — фаза);
— амплитуда;
т.е. амплитуда на выходе звена(системы) и сдвиг фазы зависят от частоты входного воздействия x(t).

Используем показательную форму записи функции единичного гармонического воздействия и отклика на это воздействие (рис. 3.1.1):

Определим связь между передаточной функцией и гармоничным воздействием, пользуясь показательной формой.
Рассмотрим звено уравнение динамики которого имеет следующий вид:

В показательной форме:

Запишем в показательной форме используя соотношения 3.1.1:

Подставим эти соотношения в (3.1.1) получим:

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

если вспомнить, что в преобразования Лапласа , то:

Получаем выражение для передаточной функции

— Амплитудно-фазовая частотная характеристика (АФЧХ)
Иногда называют частотной передаточной функцией.
Модуль АФЧХ= тождественно равен амплитуде выходного сигнала:

Сдвиг фазы выходного сигнала:

Обычно АФЧХ изображается на комплексной плоскости. Формулы (3.1.6) и (3.1.7) позволяют изобразить в полярных координатах
Так же можно изображать в традиционных декартовых координатах:

Если использовать для представления W(s) форму W(s)=K·N(s)/L(s), где L(s)- полиномы по степеням s, (причем свободные члены равны 1), а К – общий коэффициент усиления звена (системы), то

Сдвиг фазы можно определить по виду многочленов и (см. формулу (3.1.9)) т.е. как разность фаз (аргументов) числителя и знаменателя:

Постоим АФЧХ для «абстрактного» звена (системы) с передаточной функцией:

Подставляя в формулу различные значения , получаем набор векторов, на комплексной плоскости

Рассмотрим действительную и мнимую части полученных векторов Из рисунка 3.1.3 видно, что:

Амплитуда и сдвиг фазы рассчитываются для векторов, соответствующих положительным частотам и лежащих в 4 квадранте по формулам:

В общем случае для любых углов сдвига необходимо учитывать переход между квадрантами на плоскости. Тогда формула принимает вид:

где:
j = 0, 2, 3, 4. если вектор в I и IV квадрант;
j = 1, 3, 4, 4. если вектор в II и III квадранте.

Во всех технических системах отклик системы, как правило, отстает от входного воздействия, то есть сдвиг фазы всегда отрицательный. Исходя из формулы 3.1.10, степень полинома L(s) выше, чем полинома N(s). Поскольку обычно степень полинома L(s) выше, чем полинома N(s), то с увеличением частоты на входе в звено (в систему) сдвиг фазы обычно отрицателен, т.е. сигнал на выходе звена еще больше отстает по фазе от входного сигнала при увеличении частоты.
В предельном случае, если частота растет до бесконечности, мы можем вообще не получить выходного воздействия. Обычно при ω→ ∞ величина амплитуды на выходе звена стремится к 0, то есть lim A(ω→∞) = 0.

при замене на имеет зеркальное изображение.

Анализируя годографы АФЧХ при > 0 (сплошная линия на рисунке 3.1.3) и при < 0 (пунктирная линия), видим, что:
– четная функция, следовательно график симметричен относительно оси ординат, а
– нечетная функция и ее график центрально-симметричен относительно начала координат.

Кроме анализа свойств звена (системы) по годографу АФЧХ, широкое распространение получили анализ логарифмической амплитудной характеристики (ЛАХ) и фазочастотной характеристики (ФЧХ).

ЛАХ определяется как Lm(ω)=20lgA(ω).

Поскольку зачастую удобнее использовать десятичные логарифмы (lg), чем натуральные(ln), в теории управления (также и в акустике) значительно чаще используется специальная единица – децибел (1/10 часть Бела):
+1Бел – единица, характеризующая увеличение в 10 раз.
+1дБ (децибел) – соответствует увеличению в раз.

В формуле Lm(ω)=20lgA(ω) величина Lm(ω) измеряется также в децибелах. Происхождение множителя 20 таково: A(ω) – амплитуда, линейная величина, а мощность — квадратичная величина (например, напряжение в сети измеряется в Вольтах, а мощность () пропорциональна квадрату напряжения, поэтому в формуле для Lm(ω) стоит множитель 20 (чтобы привести ЛАХ (Lm(ω)) к традиционной мощностной характеристике).

Если больше на 20 дБ, то это означает, амплитуда больше амплитуды в 10 раз,

Окончательно: Lm(ω)=20lg│W(iω)│= 20lgA(ω)

Из этого следует, что +1 децибел (+1 дБ) соответствует увеличению амплитуды в раз (очень малая величина); -1 дБ – уменьшение амплитуды в раз.

Графики A(ω) и φ(ω) имеют вид:

Учитывая, что “ω” обычно изменяется на порядки и значение A(ω) – также на порядки, график Lm(ω) строится, фактически, в логарифмических координатах, т.е. Lm(ω) =Lm(lg(ω)), например:

Наклон (– 40 дБ/дек) соответствует уменьшению амплитуды в 100 раз при увеличении частоты в 10 раз.

Рассмотренные характеристики Lm(ω), то есть ЛАХ и ФЧХ имеют широкое распространение при анализе динамических свойств звена (системы), например, при анализе устойчивости САР (см. раздел “Устойчивость систем автоматического управления”).

Рисунок 3.1.10 – пример ЛАХ и ФЧХ для сложной системы

Пример 1

В качестве примера построим АФЧХ для демпфера, модель которого разобрана в этой статье. . Добавим на схему блок «Построение частотных характеристик», в качестве входа возьмем возмущающее воздействие, в качестве выхода — положение положение груза. Для наглядности иллюстрации примем в качестве выхода положение в миллиметрах (х1000), поскольку модель у нас размерная и результат получается в метрах уже достаточно маленьким примерно 0.004 метра. см. рис. 3.11

Параметры блока «Построение частотных характеристик» приведены на рисунке 3.1.12, для иллюстрации зависимости АЧХ и ЛАХ. Результат работы блока — график с выбранными параметрами — изображен на рисунке 3.1.13:

Анализ графика в линейном масштабе по ω чаще всего не очень удобен, поскольку весь график собирается в узкой области, а дальше график абсолютной амплитуды практически сливается с 0. Если мы хотим исследовать частоты хотя бы до 1000 Гц, мы увидим практически вертикальные и горизонтальные прямые. Изменения масштаба шкалы АЧХ и ω на логарифмический дает возможность лучше исследовать частотные характеристики (см. рис. 3.1.14).

На рисунке 3.1.14 представлены частотные характеристики демпфера в логарифмическом масштабе и иллюстрация соотношения между абсолютной величиной амплитуды АФЧХ и ЛАХ в децибелах.

Пример 2

Постоим частотные характеристики для чуть более сложной модели, а именно — для гидравлического демпфера, рассмотренного в предыдущей лекции.

Для начала посмотрим на модель в виде блоков.

Модель, подготовленная для анализа, представлена на рисунке 3.1.15. В отличие от исходной модели, описанной ранее, входное воздействие задается блоком «ступенька» с скачком с 0 до 1 на 10 секунде расчёта. В блоке «линейная функция» происходит пересчет сигнала «ступенька»:
0 — соответствует 200 бар в камере (конечное состояние в предыдущем примере)
1 — соответствует 400 бар в камере.
Это сделано для того, чтобы можно было подавать синусоидальный сигнал и не получать отрицательное давление в камере плунжера. Также для наглядности графика мы усиливаем выходное перемещение, переводя его из метров в миллиметры.

Частотные характеристики, получаемые в конце расчёта, приведены на рисунке 3.1.16. Видно что характеристики отличаются от простого пружинного демпфера (сравните с 3.1.14)

Блок «Построение частотных характеристик» осуществляет расчет характеристик для линеаризованной модели в окрестности заданной точки. Это означает, что частотные характеристики системы в разные моменты времени могут отличаться для нелинейных моделей. Например, в нашем случае характеристики в начале расчёта будут отличаться от характеристик, полученных в конце расчёта.

Для подробных и нелинейных моделей, блок «Построение частотных характеристик» может не работать из за наличия разрывов и нелинейностей в модели. Как например, для «точной» модели демпфера, которую мы проверяли в предыдущей статье. В этом случае возможно построить частотные характеристики непосредственно моделированием, путем подачи синусоидального сигнала с разной частотой и измерения отклика. В SimInTech для этого используется блок «Гармонический анализатор», который подключается ко входу модели и генерирует синусоидальное воздействие. В этот же блок направляется отклик системы, и производится вычисление необходимых параметров для построения различных характеристик системы, которые можно вывести на графики с помощью блока «фазовый портрет».

Модель гидравлического демпфера, собранного из библиотечных блоков SimInTech, представлена на рисунке 3.1.7

Расчеты с моделью показывают, что при сохранении общего вида графиков значения, полученные для «подробной модели», отличаются от линеаризованной модели (см. рис. 3.18 — 3.19)

Использование прямого моделирования для получения характеристик является более надежным способом и работает не только с линейными моделями, но также может быть применимо для построения характеристик некоторых реальных объектов, если их можно подключить к среде моделирования и воздействовать в реальном режиме времени. Однако затраты на вычисления значительно будут больше. Например, для получения характеристик демпфера пришлось выполнить процесс в 40 000 секунд модельного времени, на обычном компьютере это заняло порядка 35 минут. График процесса перемещения плунжера в процессе вычисления характеристик приведен на рисунке 3.1.20.

Блок «Гармонический анализатор» имеет выходы:
Re(w*t) – текущее значение действительной части амплитудно-фазовой частотной характеристики исследуемой системы;
Im(w*t) – текущее значение мнимой части амплитудно-фазовой частотной характеристики.
Это позволяет построить годограф исследуемой системы с помощью фазового портрета. (см. рис. 3.1.21)

MATLAB 8.0 (R2012b): создание, обработка и фильтрация сигналов, Signal Processing Toolbox

Окно справки по пакету расширения Signal Processing Toolbox системы MATLAB 8.0 представлено на рис. 1 на фоне рабочего окна самой системы с открытой вкладкой каталога пакетов расширения APPS. Одна из кнопок в панели инструментов Signal Analysis дает доступ к браузеру сигналов, фильтров и спектров, показанному в правой части окна справки. В левой части этого окна указаны наименования разделов пакета Signal Processing Toolbox.

Рис. 1. Окно справки по пакету расширения Signal Processing Toolbox

Как видно на рис. 1, пакет Signal Processing Toolbox состоит из следующих разделов:

• Waveforms — создание сигналов с различной формой и разными законами модуляции;
• Convolution and Correlation — свертка и корреляция сигналов ;
• Transform — преобразование сигналов ;
• Analog and Digital Filters — аналоговые и цифровые фильтры ;
• Spectral Analysis — спектральный анализ.

Создание сигналов

Многие сигналы представлены как функции времени s(t), параметры которой можно изменять с помощью модуляции того или иного вида. Модуляцией называют процесс изменения какого-либо параметра (амплитуды, частоты, фазы и т. д.) по определенному закону, в результате чего сигнал становится переносчиком информации.

MATLAB 8.0 со своими встроенными средствами позволяет создавать множество сигналов. Например, простейшим является синусоидальный сигнал:

s = A sin(2π ft +j),

где A — амплитуда; f — частота; j — фаза сигнала. Задав эти параметры ( t — как вектор отсчетов сигнала, число элементов которого определяет число отсчетов сигнала):

можно легко построить график сигнала s ( t ):

Он показан на рис. 2 в графическом окне системы MATLAB. Такое построение осуществляется в командной строке, на что указывает приглашение к вводу >>. В командном окне также показано окно с краткими данными о системе MATLAB 8.0 (R2012b). Согласно этим данным система выпущена на рынок в 2012 году.

Рис. 2. Окно командного режима

В дальнейшем, учитывая множество сигналов, мы будем представлять их по четыре в каждом графическом окне, используя для этого функцию его разбиения на четыре под-окна subplot.

Рис. 3. Окно с подокнами меандра, апериодического треугольного импульса, пилообразного периодического импульса и периодического треугольного импульса

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

Рис. 4. Временные диаграммы четырех сигналов, относящихся к функциям Гаусса

Важное значение имеет функция sinc (или sin(π t )/π t при t ≠ 0 и 1 при t = 0). Функция sinc( t ) представляет обратное преобразование Фурье для прямоугольного импульса с высотой 1 и шириной 2π:

Кроме того, эту функцию можно использовать как базисную для восстановления любого сигнала g ( t ) по его отсчетам, если спектр сигнала ограничен условием –p < w < p:

Рис. 5. Получение непрерывной кривой, проходящей через точки (отсчеты) произвольного сигнала

Это положение, вытекающее из известной теоремы Котельникова, иллюстрирует приведенный ниже пример для десяти случайных точек сигнала (рис. 5):

Подобный метод восстановления аналогового сигнала из его цифровых отсчетов ныне применяется во всех цифровых осциллографах.

Рис. 6. Сложные модулированные сигналы

Сложные модулированные сигналы показаны на рис. 6, здесь chirp — период частотно-модулированного сигнала, sawtooth — сигнал с частотной модуляцией по треугольному закону, pulstrain — скачок с линейным спадом, pulstran gauspuls — последовательность импульсов Гаусса с убывающей амплитудой. Для создания этих четырех модулированных сигналов служит программа:

формирует выборку (дискретные значения) косинусоидального сигнала с частотой от f 0 в начальный момент времени t до f 1 в конечный момент времени t 1. Звук такого сигнала напоминает визг, откуда и его название (chirp). По умолчанию t = 0, f 0 = 0 и f 1 = 100. Необязательный параметр phi (по умолчанию 0) задает начальную фазу сигнала. Другой необязательный параметр — method — задает закон изменения частоты: linear — линейный закон (по умолчанию), quadratic — квадратичный и logarithmic — логарифмический.

Функция strips позволяет детально (по частям с длительностью 0,25 с) рассмотреть первые два сложных сигнала, например частотно-модулированного от функции vco. В разделе модуляция/демодуляция есть пример с GUI (рис. 7), иллюстрирующий такие сигналы с различными видами модуляции.

Рис. 7. Визуализация модулированных сигналов

Свертка и корреляция сигналов

Пусть имеется две последовательности, представленные векторами a и b . Сверткой называют одномерный массив, вычисляемый следующим образом:

При записи этого выражения учтено, что нумерация индексов массивов в MATLAB идет с единицы (1). Свертка реализуется функцией conv(a,b):

Операцию свертки часто используют для вычисления сигнала на выходе линейной системы y по сигналу на входе x при известной импульсной характеристике системы h :

Эту операцию можно использовать для осуществления простейшей фильтрации сигнала:

Для двух векторов x и y с длиной m и n определена операция свертки:

Рис. 8. Линейная и циклическая свертка

Обратная свертке функция — это [q,r] = deconv(z,x). Она фактически определяет импульсную характеристику фильтра (рис. 8):

Рис. 9. Очистка от шума линейно-нарастающего сигнала

Свертку часто применяют для очистки сигналов от шума (рис. 9):

Для двумерных массивов также существует функция свертки: Z = conv2(X,Y) и Z = conv2(X,Y,‘option’). Возможна и многомерная свертка — функция convn. Новая функция mscohere строит график зависимости квадрата модуля функции когерентности от частоты (рис. 10):

Рис. 10. Зависимость квадрата модуля функции когерентности от частоты

Следующая программа строит отсчеты двух сигналов с задержкой на три отсчета — треугольного сигнала и шума (рис. 11):

Рис. 11. Два коррелированных сигнала — треугольный и шума

Кросс-корреляцию этих сигналов обеспечивает программа (рис. 12):

Рис. 12. Кросс-корреляция сигналов

Преобразование сигналов

В MATLAB всегда много внимания уделялось различным преобразованиям сигнала [2–5], в частности спектральным. Спектр дискретного сигнала является периодическим, и прямое дискретное преобразование Фурье (ДПФ или Discrete Fourier Transform, DFT) определяется выражением:

Для предотвращения растекания (размазывания) спектра дискретных сигналов часто используются окна. Для этого достаточно в формуле прямого ДПФ под знаком суммы ввести еще один множитель — W ( k ). Соответственно, обратное дискретное преобразование Фурье задается выражением:

ДПФ легко обеспечивает восстановление непрерывных периодических сигналов с ограниченным спектром. Для этого нужно номер отсчета k поменять на нормированное время t / T . Тогда формула восстановления при четном числе отсчетов будет иметь вид:

Для получения полосы частот сигнала от 0 до π/ T приходится смещать нумерацию отсчетов. При нечетном числе отсчетов суммирование ведется при n , меняющемся от –( N –1)/2 до ( N –1)/2. Коэффициенты X ·( n ) с отрицательными номерами вычисляют из соотношения симметрии.

Частотным спектром случайного процесса является преобразование Фурье от корреляционной функции случайного процесса R_x :

В радиоэлектронике особый интерес представляет спектральная оценка сильно зашумленных сигналов. Для таких сигналов применяются два подхода: непараметрический — использующий только информацию, извлеченную из сигнала (реализован в методах периодограмм и Уэлча), и параметрический — предполагающий наличие некоторой статистической модели сигналов, параметры которой подлежат определению. Реализовано восемь классов алгоритмовспектрального анализа : Periodogram, Welch, MTM (Thomson multitaper method), Burg, Covariance, Modified Covariance, Yule-Walker, MUSIC (Multiple Signal Classification) и Eigenvector. Их подробное описание дано в [4–6].

Периодограммы, спектрограммы и их применение

Обычный спектр строится методом быстрого преобразования Фурье (БПФ, FFT) часто с применением временного окна, предотвращающего разрывы сигнала на концах интервала анализа спектра. Например, так реализованы периодограммы. Для построения спектрограмм используется разбивка интервала анализа на короткие окна. Короткое окно пробегает общий интервал анализа, и в каждом частичном интервале строится свой спектр. Их наложение дает спектрограмму, определенную в пространстве «уровень – частота – время» (на плоскости уровень представляется цветом), тогда как обычный спектр определен в плоскости «уровень – частота».

Рис. 13. Построение спектрограмм сигналов, модулированных по различным законам: логарифмическому, линейному, треугольному и синусоидальному

На рис. 13 показано построение четырех спектрограмм для сигналов с различными законами частотной модуляции. Нетрудно увидеть, что во всех случаях закон модуляции отчетливо распознается и позволяет судить об области частот, в которой действует сигнал (chirp или vco). Этот рисунок строит следующая программа:

Рис. 14. Сравнение периодограммы и спектрограммы сигнала sinc(t)

Из сказанного может сложиться неверное представление о явных преимуществах спектрограмм по сравнению с периодограммами (функция psd). То, что это далеко не так, показывает программа для сигнала sinc (рис. 14):

Эта программа строит периодограмму и спектрограмму функции sinc(t). Теперь беспомощной оказывается спектрограмма, по которой ничего нельзя сказать о спектре сигнала и области его частот. А периодограмма более информативна: она указывает на вид спектра и область занимаемых им частот. В частности, хорошо видно постоянство спектра в начальной области частот, присущее этой функции. Правда, вид спектрограммы сильно зависит от типа короткого окна: в данном случае задано окно Блэкмана (Blackman). Вид периодограммы также зависит от выбора окна, но глобального.

Чем сложнее сигнал, тем более детальной и эффектной оказывается спектрограмма. На рис. 15 показан пример модуляции/демодуляции с построением спектрограммы сложного звукового сигнала. Кстати, оригинальный сигнал и сигнал, прошедший модуляцию/демодуляцию, можно воспроизвести на компьютере, оборудованном звуковой картой и акустической системой. Это можно сделать при различных видах модуляции.

Рис. 15. Спектрограмма сложного акустического сигнала с амплитудной модуляцией

Похожие на спектрограммы картинки дают вейвлетограммы и скайлеграммы, получаемые при вейвлет-анализе сигналов [7]. Порою они более информативны. Но вейвлеты в Signal Ptocessing Toolbox не реализованы: они описаны и используются в отдельном пакете расширения Wavelet Toolbox.

Для выполнения дискретного преобразования Фурье различными методами служит GUI-окно Discrete Fourier Transform (рис. 16).

Рис. 16. Спектр пилообразного сигнала в окне дискретного Фурье-преобразования

Оно позволяет задать один из трех видов сигнала (синусоида, меандр и пилообразный) и его периодограмму (спектр) при одном из шести видов окон. Как сигнал, так и окно можно загружать извне из файла или рабочего пространства MATLAB. Нетрудно убедиться в большом влиянии на вид спектра выбранного окна.

Оконные функции и браузер окон

Учитывая важную роль окон, в Signal Processing Toolbox входит 21 N ‑точечное окно ( N — целое число), называемое по фамилии предложивших окно ученых, например Hamming, Blackman, Bartlett, Chebyshev, Taylor, Kaiser и т. д. Для просмотра временных и амплитудно-частотных характеристик всех окон есть соответствующие функции, но удобно пользоваться GUI-браузером окон wvtool.

позволяет строить четыре типа 64‑точечных окон — прямоугольное, Хемминга, Ханна и Гаусса. Можно задать окна и отдельными командами.

Рис. 17. Сравнение в браузере трех окон Кайзера с параметром бетта 1,5, Блэкмана и Блэкмана-Харисса

Это, а также сравнительное построение трех других типов окон обеспечивает следующая программа (рис. 17):

Рис. 18. Конструктор окон

Существует также конструктор-анализатор окон с GUI-интерфейсом, окно которого (рис. 18) открывается командой:

По умолчанию в нем открывается окно объекта sigwin, но можно открыть и другие окна (в том числе окно пользователя) из списка Current Window Information. Открытые окна появляются в третьем нижнем окне. Окна можно скопировать, добавить в список, установить в рабочее пространство MATLAB или стереть.

Изменение числа отсчетов и интерполяция сигналов

Изменение числа отсчетов широко используется в технике цифровой обработки сигналов. Для уменьшения числа отсчетов применяются операция децимации и функция decimate (рис. 19):

Рис. 19. Уменьшение частоты дискретизации сигнала (децимация)

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

Рис. 20. Интерполяция сигнала и увеличение частоты его дискретизации

Аналоговые и цифровые фильтры

Фильтры имеют особое значение при обработке сигналов. С их помощью осуществляется очистка сигналов от шума или реализуются избирательные свойства систем. Будем считать, что читатель знаком с теорией фильтров.

Фильтрующие цепи обычно задаются своей операторной передаточной характеристикой:

h ( s ) = a ( s )/ b ( s ).

Имея векторы коэффициентов полиномов a ( s ) и b ( s ), с помощью функции freqs можно построить АЧХ и ФЧХ фильтрующей цепи в логарифмическом масштабе (рис. 21):

Рис. 21. Построение логарифмической АЧХ и ФЧХ системы по ее операторной передаточной характеристике

В Signal Processing Toolboox входит множество функций по расчету и проектированию различных фильтров — нижних, верхних частот и полосовых. Порою для получения важных характеристик фильтров достаточно задать нужную строку программного кода. Например, построение характеристик аналогового фильтра Бесселя (рис. 22) реализуется следующим программным фрагментом:

Для просмотра практически всех характеристик фильтров можно использовать визуализатор фильтров fvtool. Приведенная ниже программа дает характеристики фильтра Бесселя нижних частот (рис. 23):

Рис. 23. Характеристики аналогового эллиптического фильтра нижних частот

Открытое меню Analysis дает представление обо всех доступных видах анализа фильтров. Среди них построение АЧХ и ФЧХ фильтра, импульсная и переходная характеристики, групповая задержка и т. д.

Конструирование двух цифровых фильтров НЧ (FIR и Баттерворта) со сравнительным построением их характеристик в окне fvtool (рис. 24) обеспечивает следующая программа:

Рис. 24. АЧХ двух цифровых фильтров для сравнения

Программа строит и другие характеристики фильтров, например переходные и импульсные. Напомним, что переходная характеристика является реакцией системы (фильтра) на единичный скачок, а импульсная — на импульс единичной площади с длительностью, стремящейся к нулю (рис. 25).

Рис. 25. Импульсные характеристики двух цифровых фильтров

Интерактивный конструктор-анализатор фильтров fdatool

В Signal Proccesing Toolbox входит интерактивный конструктор-анализатор фильтров с GUI-интерфейсом. Он позволяет без какого-либо программирования анализировать и проектировать семь основных типов фильтров (таблица).