Тема: «Спектральный анализ сигналов»
Курсовая работа выполняется самостоятельно на заключительном этапе изучения курса и ставит целью освоение возможностей визуального программирования и прикладных пакетов для вычислений и визуализации данных при решении учебных и профессиональных задач. Для выполнения курсовой работы требуется:
· изучить язык программирования (BASIC или ПАСКАЛЬ), типы данных, основные операторы, функции и процедуры, механизм передачи параметров в них;
· изучить основные алгоритмы вычислительной математики: решение уравнений (поиск корня уравнения), численное интегрирование, аппроксимация данных по методу наименьших квадратов, численное решение обыкновенных дифференциальных уравнений, методы гармонического анализа периодических функций;
· освоить приемы работы в визуальных средах программирования (Visual Basic, Delphi),использование объектов, их свойств и событий;
· освоить приемы работы в математических (Maple) и инженерных (MathCAD, MATLAB) пакетах для обработки и визуализации данных и результатов вычислений;
· проделать вычисления демонстрационного варианта программы на языке системы Maple для сигнала своего варианта.
Постановка задачи:
Написать программу на языке программирования (Бейсик, Паскаль) для решения следующей задачи (10 вариантов задания). Построить блок-схемы задачи и вспомогательных частей алгоритма. Оформление графиков и таблиц выполнять средствами математических и инженерных пакетов.
Задача:
Рассматривается функция F(t), представляющая собой периодический (период 2π) сигнал единичной амплитуды длительности T. Значение величины T есть наименьший положительный корень полинома (таблица 1- по вариантам), который вычисляется любым из известных методов нахождения корней уравнений.
1. Программу, реализующую вычисление корня Т, написать на языке программирования Pascal ABC. Для построенного сигнала вычислить по формулам Бесселя коэффициенты конечной суммы Фурье и записать их в файл.
2. С помощью языка программирования математического пакета Maple вычислить для аналогового сигнала своего варианта коэффициенты конечной суммы Фурье (аппроксимирующего тригонометрического полинома) в аналитической форме, построить графики исходного сигнала, тригонометрического полинома для нескольких степеней, а также графики спектров амплитуд дискретного (с коэффициентами из файла, вычисленными на Паскале) и аналогового сигнала.
Оформление:
v Формат А4.
v Титул
v Постановка задачи
v Алгоритмы решения вспомогательных задач
v Блок-схемы
v Общая структура программы на языке программирования
v Результаты расчетов, графики
v Литература
Таблица 1. Варианты полиномов
| № варианта | Полином |
| x^5-8*x-1 | |
| x^5-3*x-3 | |
| 2*x^4-x^3-8 | |
| x^6-4*x^4-2 | |
| x^6-4*x^4-3 | |
| x^6-3*x^3-2 | |
| x^5-7*x-14 | |
| x^5-x^3-1 | |
| x^5-2*x^3-4 | |
| x^6-3*x^4-5 |
Решение задачи средствами я.п. Pascal.
uses CRT;
const
epsilon=0.00000001;
Npoint=31;{число точек дискретизации}
nkoef=50;
h=2*Pi/Npoint;{шаг дискретизации}
type
fourier=array[0..nkoef] of real;
VAR
a,b,ya,yb,x,y:real;
x0,delta,x1:real;
ak,bk:fourier; {коэффициенты}
T:real; {корень полинома }
function f(t:real):real;{полином}
begin
f:=t*t*t*t*t-2*t*t*t-4; {здесь нужна формула своего варианта полинома}
{f:=t*t*t*t*t-t*t*t-1;}
end;
function f1(t:real):real;{производная полинома}
begin
f1:=5*t*t*t*t-6*t*t;
end;
function signal(x,T:real):real;
var
y:real;
begin
y:=0;
if x<=T then y:=1;
signal:=y;
end;
procedure Bessel(m,Nt:integer;var A,B:fourier);
var
p,q:real;
k,j:integer;
begin
for k:=0 to m do
begin
p:=0;q:=0;
for j:=0 to Nt-1 do
begin
x:=j*h; {h-- шаг дискретизации}
y:=signal(x,T); {T -глобальная переменная}
p:=p+y*cos(k*x);
q:=q+y*sin(k*x);
end;
A[k]:=2*p/Nt;
B[k]:=2*q/Nt;
writeln(k:4,A[k]:8:4,B[k]:12:4);
end;
end;
procedure zapis(m:integer;A,B:fourier);
var
k:integer;
fn:string;
ft:text;
begin
fn:='D:\коэфф.txt';
assign(ft,fn);
rewrite(ft);
for k:=0 to m do
writeln(ft,k:6,A[k]:10:4,B[k]:10:4);
close(ft);
end;
{ ===========================================}
begin
clrscr;
{----- локализация корня------}
repeat
writeln('Задайте a,b');
readln(a,b);
ya:=f(a);
yb:=f(b);
writeln(a:6:2,b:5:1,ya:12:8,' ',yb);
until ya*yb<0;
{==========метод касательных=====================}
x0:=b;
delta:=1;
repeat
x1:=x0-f(x0)/f1(x0);
delta:=abs(x0-x1);
x0:=x1;
writeln (x0);
until delta<epsilon;
{==============================================}
writeln('корень=',x0,y:18:8);
readln;
T:=x0;
{вычислен корень Т}
{-------- коэффициенты Фурье по формулам Бесселя---- }
Bessel(nkoef,Npoint,ak,bk);
readln;
zapis(nkoef,ak,bk);
end.
Результаты вычислений
0 0.5806 0.0000
1 0.3474 0.3655
2 -0.0157 0.3099
3 -0.0645 0.0554
4 0.0790 0.0080
5 0.0804 0.1039
6 -0.0124 0.0809
7 0.0082 -0.0057
8 0.0740 0.0152
9 0.0405 0.0649
10 -0.0057 0.0220
11 0.0359 -0.0199
12 0.0656 0.0206
13 0.0191 0.0390
14 0.0046 -0.0123
15 0.0533 -0.0229
16 0.0533 0.0229
17 0.0046 0.0123
18 0.0191 -0.0390
19 0.0656 -0.0206
20 0.0359 0.0199
21 -0.0057 -0.0220
22 0.0405 -0.0649
23 0.0740 -0.0152
24 0.0082 0.0057
25 -0.0124 -0.0809
26 0.0804 -0.1039
27 0.0790 -0.0080
28 -0.0645 -0.0554
29 -0.0157 -0.3099
30 0.3474 -0.3655
31 0.5806 0.0000
32 0.3474 0.3655
33 -0.0157 0.3099
34 -0.0645 0.0554
35 0.0790 0.0080
36 0.0804 0.1039
37 -0.0124 0.0809
38 0.0082 -0.0057
39 0.0740 0.0152
40 0.0405 0.0649
41 -0.0057 0.0220
42 0.0359 -0.0199
43 0.0656 0.0206
44 0.0191 0.0390
45 0.0046 -0.0123
46 0.0533 -0.0229
47 0.0533 0.0229
48 0.0046 0.0123
49 0.0191 -0.0390
50 0.0656 -0.0206
Решение задачи средствами Maple
Программа построения АЧХ для периодической функции, заданной аналитически на периоде
> restart;with(plots):;with(plottools):;
> N:=51;
> P:=x^5-2*x^3-4; #полином
> T:=fsolve(P,x=0..3);
N - число гармоник
> plot(P,x=T-0.5..T+0.5,thickness=2,color=blue):;
> 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):;
>
Аппроксимация функции конечной суммой ряда Фурье есть тригонометрический полином
степени 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;
>
>
>
>
> 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);
num[0]:=0;
> 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,symbolsize=20,color=blue,title=`спектр амплитуд`):
> display(Ar,Rsym);
Чтение коэффициентов Фурье, вычисленных программой на Паскале
> fn:=`d:\\коэфф.txt`;
> L:=readdata(fn,3):;
> with(linalg):
> 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,Ar,Rsym);
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>