5. Спектральный анализ программы средствами математического пакета maple

> restart;

> with(linalg):with(plots):

pp:=(x,y)->[x,y];

> f:= proc(t) local z ;z:=piecewise(t<0,0,t<=1,1,0); evalf(z);end;

> plot(f(t),t=-1..2,thickness=2,color=brown);

> T:=2*Pi;#период сигнала

> tau:=2;#длительность сигнала

Для полученной периодической функции вычисляем коэффициенты тригонометрического ряда Фурье

в соответствии с формулами

Коэффициенты Фурье вычисляются в аналитическом виде процедурой Fourier(T,F), где параметрами процедуры являются величина периода T и имя процедуры, вычисляющей значения сигнала на отрезке периодичности. В нашем случае это процедура Period, для использования, которой следует согласовать параметры обеих процедур. Процедура Fourier вычисляет коэффициенты Фурье и передает их как глобальные для процедуры объекты. Интегрируемая функция имеет один параметр - переменную интегрирования. При этом функция должна быть определена на отрезке периодичности, поэтому следует задать процедуру, в которой не должно быть других параметров:

> signal:=proc(x)global tau, T;

f(x/tau);

end;

> plot(signal(x),x=0..T,thickness=2,color=blue);

Таким образом, фактическим параметром, соответствующим названию функции при вызове процедуры Fourier должна быть процедура F.

> Fourier:=proc(T,f) global a0,A,B;

a0:=2/T*evalf(Int(f(x),x=0..T)):;

A:=value(2/T*int(f(x)*cos(2*Pi*k*x/T),x=0..T));

B:=value(2/T*int(f(x)*sin(2*Pi*k*x/T),x=0..T));

end proc:

Вычислим выражения для коэффициентов и выведем их значения:

> Fourier(T,signal):

> a0;evalf(A,4);evalf(B,4);

Значения коэффициентов и можно вычислить подстановкой (процедура subs)

натуральных чисел в полученные выражения

Зафиксируем число коэффициентов (степень )тригонометрического полинома

Фурье и построим совместные графики исходного сигнала и приближающих его

полиномов для нескольких значений степеней. Глобальный параметр k, как и выражения для коэффициентов

Фурье определяются процедурой Fourier

> Trig:=proc(t,n) local z;global k,T,a0,A,B;

z:=a0/2+sum(A*cos(k*2*Pi*t/T)+B*sin(k*2*Pi*t/T),k=1..n);

end:

Построим график триг. полинома пятой степени

> plot(Trig(x,5),x=-T..2*T);

> grafik:=plot(signal(x),x=-0.2..T+0.2,thickness=3,color=blue,discont=true):

>

Построим совместные графики сигнала и полиномов нескольких степеней

> CL:=brown,red,gold;

> RisTrig:=seq(plot(Trig(t,2*j),t=-0.1..2*T,color=CL[j]),j=1..3):

> display(RisTrig,grafik);;

Часть 2

> with(plottools):

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

> Fourier_T:=proc(F,T0,TN,k::evaln) local T;

global A0,Ak,Bk;

T:=TN-T0;

A0:=2/T*Int(F(x),x=T0..TN);

Ak:=2/T*int(F(x)*cos(k*x*2*Pi/T),x=T0..TN):

Bk:=2/T*int(F(x)*sin(k*x*2*Pi/T),x=T0..TN):

end proc:

Процедура Trig_polynom подстановкой номеров k в формулы коэффициентов вычисляет эти коэффициенты и сохраняет во вновь созданных массивах a,b [0..N]. Заодно строится график тригонометрического полинома степени N.

> Trig_polynom:=proc(N,T,a0,ak,bk,k::evaln) local n,Pol,A0,A,B;

global a,b,RisTrig;

a:=array(0..N);b:=array(0..N);

A0:=evalf(a0);a[0]:=A0;b[0]:=0;

A:=seq(evalf(subs(k=n,ak)),n=1..N);

B:=seq(evalf(subs(k=n,bk)),n=1..N);

for n from 1 to N do

a[n]:=A[n];b[n]:=B[n];

end do;

Pol:=A0/2+sum(A[k]*cos(2*Pi*k*x/T)+B[k]*sin(2*Pi*k*x/T),k=1..N):

RisTrig:=plot(Pol,x=-T/2..3*T/2,color=blue,thickness=2):

RETURN(Pol);

end proc:

Процедура ARR строит спектр амплитуд линейчатой гистограммой.

> ARR:=proc(n::integer,c) local L,H,ma,mi,k::integer,

Sim::array(0..n);

ma:=c[0];mi:=c[0];

L:=line([0,c[0]],[n,c[n]],thickness=2,color=red);

for k from 1 to n do

if c[k]>ma then ma:=c[k];end if;

if c[k]<mi then mi:=c[k];end if;

end do;

H:=ma-mi;

if H=0 then RETURN(L) end;

for k from 0 to n do

if abs(c[k])<H/1000 then

Sim[k]:=ellipse([k,c[k]],0.2,0.01*H,color=blue);

else

Sim[k]:=plottools[arrow]([k,0],[k,c[k]],0.2,0.2,0,color=black);

end if;

end do;

convert(Sim,list);

end:

Процедура Spectr вычисляет модули и аргументы комплексной (синусной) формы тригонометрического полинома.

График спектра фаз - структура Risphi.

> Spectr:=proc(n,a,b,c,Risphi) local k,R,phi;

for k from 0 to n do

c[k]:=evalf(abs(I*a[k]+b[k])):

phi:=evalf(argument(I*a[k]+b[k]));

R[k]:=[eval(k),eval(phi)];

end:;

Risphi:=plot(convert(R,list)):

end:

Определим прямоугольный импульс f(t) единичной длительности.

Построим сигнал (F_for_all) длительности и для него выполним вычисления спектральных характеристик.

Здесь и дальше команды не будут зависеть от варианта сигнала, поскольку название обрабатываемого сигнала теперь у всех одинаково (F_for_all)

Положим величину периода T=2 (для всех вариантов). Аналитические выражения коэффициентов вычисляются в процедуре Fourier_T.

> tau:=Pi;# длительность сигнала

> T:=2*Pi;# величина периода

> F_for_all:=proc(t) global tau;f(t/tau);end proc:;

> Ris1:=plot(F_for_all(t),t=0..T,color=brown,thickness=2):display(Ris1);

Получим формулы для коэффициентов Фурье, зависящие от номера коэффициента k.

> Fourier_T(F_for_all,0,T,k):

> a0:=evalf(A0);a_k=Ak;b_k=Bk;

Значения величин С[k] в аналитическом виде можно вычислить (для любого натурального k) по формуле:

> C[k]:=simplify(sqrt(Ak^2+Bk^2));

Для построения спектра амплитуд зафиксируем число коэффициентов и с помощью процедуры Trig_polynom вычислим коэффициенты a[k], b[k] тригонометрического полинома степени N, а также построим его график на отрезке, большем длины периода. Массивы коэффициентов создаются при работе процедуры и являются глобальными, то есть доступными вне процедуры.

> N:=19;

> Trig_polynom(N,T,A0,Ak,Bk,k):

> display(RisTrig,Ris1);

Амплитудный спектр (коэффициенты С[k]) вычисляется процедурой Spectr, которая как побочный эффект строит графическую структуру Risphi1, являющуюся графиком спектра фаз. График спектра амплитуд строит процедура ARR.

> Spectr(N,a,b,c,'Risphi1'):;

> display(ARR(N,c));

> RisAMP:=display(ARR(N,c)):

> display(Risphi1);