Спектральный анализ сигналов
Ниже приведены команды системы математических вычислений MAPLE с комментариями и минимальными пояснениями используемых алгоритмов.
Команда restart нужна для очистки памяти системы при повторном запуске программы и при начале следующего сеанса. Следующие команды подключают пакеты линейной алгебры и пакеты графики, соответственно. Точка с запятой, завершающая команду, служит признаком отображения результата выполнения команды. Двоеточие подавляет вывод результата (в процедуре --- это последний выполненный оператор).
> restart;
> with(linalg):with(plots):with(plottools):
Полином
Здесь следует задать свой вариант полинома, наименьший положительный корень которого используется как длительность T сигнала значение которого вычисляется командой fsolve (решить уравнение p(x)=0 относительно переменной x в плавающем формате с начальным приближением из отрезка [0;3] ). Если не указать вариант в плавающей форме (команда solve), то система будет пытаться решить уравнение в аналитическом виде. Команда plot строит график на интервале [T-0.5, T+0.5].
> p(x):=x^5-8*x-1;
> Koeff:=fsolve(p(x)=0,x,0..3);
> tau:=Koeff:
> plot(p(x),x=Koeff-0.5..Koeff+0.5,thickness=2,color=black);
Далее строится сигнал длительности Т единичной амплитуды.
> f:=proc(t) local z;
z:=piecewise(t<0,0,t<T,1,0);end:
> plot(f(x),x=-1..2*Pi,color=blue,thickness=2,discont=true);
Коэффициенты ряда Фурье вычисляются в аналитической форме
(если функция F(t) допускает вычисление первообразной) по формулам
В противном случае интегралы вычисляются встроенными алгоритмами численного интегрирования. Интегрирование выполняется по любому отрезку длины периода:
>
> a0:=1/Pi*Int(f(t),t=-Pi..Pi):value(%);
>
> ak:=1/Pi*Int(f(t)*cos(k*t),t=-Pi..Pi):value(%);
> bk:=1/Pi*Int(f(t)*sin(k*t),t=-Pi..Pi):;value(%);
Вычисляем последовательность коэффициентов Фурье для построения спектра амплитуд аналогового сигнала.
> A:=seq(evalf(subs(k=n,ak),5),n=1..N);
> B:=seq(evalf(subs(k=n,bk),5),n=1..N);
> C:=seq(evalf(sqrt(A[n]^2+B[n]^2),3),n=1..N);
Аппроксимация функции конечной суммой ряда Фурье есть тригонометрический полином степени n.
> Trig:=proc(t,n) local z,k;global a0,A,B;
z:=a0/2+sum(A[k]*cos(k*t)+B[k]*sin(k*t),k=1..n);
end;
Процедура Trig вычисляет значения тригонометрического полинома в точке t, что дает возможность построить график полинома степени N.
> plot(Trig(x,10),x=-2*Pi..2*Pi,numpoints=1000);
Меняя число коэффициентов, строим совместный график сигнала и нескольких тригонометрических полиномов.
Совместный график:
> ris1:=plot(Trig(x,20),x=-2*Pi..2*Pi,numpoints=1000):
> ris2:=plot(f(x),x=-Pi..Pi,thickness=3,color=blue):
> display(ris1,ris2);
> col:=[brown,black,green];
> Ris:=seq(plot(Trig(x,3*n),x=-2*Pi..2*Pi,numpoints=100,color=col[n]),n=1..3):
> display(Ris,ris2);
График спектра амплитуд аналогового сигнала
> k:='k';c:=array(0..N);num:=array(0..N);
> c[0]:=evalf(abs(a0),3);
> for k from 1 to N do
c[k]:=evalf(abs(A[k]+I*B[k]),3):
num[k]:=k;
#print(k,num[k],c[k]);
end:;
>
> Risc:=zip((x,y)->[x,y],num,c):
> arr:=array(0..N);
> for i from 0 to N do
arr[i]:=arrow([i,0],[i,c[i]],0.2,0.2,0,color=black);
end:;
> Ar:=convert(arr,list):;
> Rsym:=plot(Risc,style=point,symbol=circle,thickness=2,color=black,title=`спектр амплитуд`):
> display(Ar,Rsym);
Дискретизация задачи
Для построения спектральных характеристик дискретного сигнала в Maple используем результаты расчетов, выполненные на Паскале (временный файл с коэффициентами Фурье).
> fn:=`D:\\koeff`;
> L:=readdata(fn,3):;
> with(linalg):
Warning, the protected names norm and trace have been redefined and unprotected
> Nkoef:=vectdim(L);
> j:='j';
> for j from 1 to Nkoef do
c[j-1]:=evalf(sqrt(L[j,2]^2+L[j,3]^2),4);
#print(j,c[j-1]);
end:
Строим спектр амплитуд дискретного сигнала и сравниваем графики этих спектров.
> for i from 0 to Nkoef do
arr[i]:=arrow([i,0],[i,c[i]],0.2,0.2,0,color=blue);
end:;
> Ardiskr:=convert(arr,list):;
> Rlin:=plot(Risc,thickness=2,color=brown,title=`сравнение спектров амплитуд`):
> display(Ardiskr,Rlin);
>
Осмысление результатов.
На этом этапе полезно выполнить расчеты для разных значений параметров Nt. (Nt < 100) – число точек дискретизации и -- число коэффициентов тригонометрического полинома. Обратите внимание на эффект Гиббса вблизи точек разрыва аппроксимируемой функции Это на совместном графике аналогового сигнала и тригонометрического полинома.
Литература
Основная:
1 Кондратьев В.П. Вычислительная математика. Учебное пособие. Екатеринбург: УрТИСИ, 2010.
2 Рыжиков Ю.И. Вычислительные методы. СПб.: БХВ-Петербург, 2007.
3 Амосов А.А., Дубинский Ю.А., Копченова Н.В. Вычислительные методы для инженеров.-М.: Высшая школа, 1994
4 Бахвалов Н.С. Численные методы. -М.: Наука, 1973
5 Фаронов В.В. DELPHI. Программирование на языке высокого уровня. Учебник для вузов.—СПб.:Питер, 2005, 640с.
6 Немнюгин С. Turbo Pascal. Программирование на языке высокого уровня. Учебник для вузов. 2-е изд.—СПб.:Питер, 2005, 544с.
7 Немнюгин С. Turbo Pascal. Практикум. Учебное пособие. 2-е изд.—СПб.:Питер, 2004, 272с.
8 Кондратьев В.П. Языки программирования. Система Maple.ч. I. Основы работы в системе. Ч. III/. Язык программирования системы. Учебное пособие. Екатеринбург: УрТИСИ ГОУ ВПО «СибГУТИ», 2006.
Дополнительная:
9 Д. Мак-Кракен, У. Дорн. Численные методы и программирование на Фортране. – М.: Мир, 1977.
10 Хемминг Р.В. Численные методы. – М.: Наука, 1972.
11 Дьяконов В.,П. Maple 7: учебный курс. – СПб.: Питер, 2002.
12 Сдвижков О.А. Математика на компьютере: Maple 8. – М.: СОЛОН-Пресс, 2003.
13 Говорухин В.Н., Цибулин В.Г. Компьютер в математическом исследовании. Учебный курс. СПб: Питер, 2001
Говорухин В.Н., Цибулин В.Г. Компьютер в математическом исследовании. Учебный курс. 2001