Mathcad - ЛР4-2

Найдем коэффициенты с, решив СЛУ А*с=B

X := lsolve(A,B)	æ-0.380449692512293 ö
	ç	0.86314784758605	÷
	ç		÷

X= ç-0.988881427117785 ÷ çç 3.68361142780256 ÷÷ è -8.59815999932724 ø

			æ	0	ö
			ç	-0.380449692512293	÷
C(X ,n) :=	for i Î 0.. n - 2		ç		÷
	Ci+1	¬ Xi	ç	0.86314784758605	÷
			c = ç	-0.988881427117785	÷
	Cn ¬	c := C(X ,n)
		0	ç	3.68361142780256	÷
	C0 ¬ 0		ç		÷
			ç	-8.59815999932724	÷
			ç		÷
	C		è	0	ø

А также b&d:

b(h,g, c,n) := for i Î 0.. n - 1

b ¬ g - 1 × h × éc + + (2 × c) ù i i 3 i ë i 1 iû

b

b := b(h, g,c, n)

d := d(h,c, n)

æ0.668165641707643 ö çç-3.13633128341529 ÷÷

ç-1.96643804814484 ÷

b= ç-4.12317048118012 ÷ çç 14.4103123812407 ÷÷ è 16.8629866653211 ø

d(h,c, n) :=	for	i Î 0.. n - 1
	d	¬	ci+1 - ci
	i		3 × hi
	d

æ-0.012681656417076 ö çç 0.103633128341529 ÷÷ ç-0.077167886445993 ÷

d = ç 0.222499659758112 ÷ çç -4.09392380904327 ÷÷ è 0.95535111103636 ø

Составим общий вид многочлена на заданном отрезке t:

SP(z) := t ¬ n while z < xt

t ¬ t - 1

at + bt × (z - xt) + ct × (z - xt)2 + dt × (z - xt)3

Найдем пропущенное значение для х= -5:

SP(−5) = −28.0976277975134

Выведем на один чертеж график кубического сплайна и точечный график функции.

i := 0.. n

		15
20	10	1.25	0	10
SP(z)
yi		17.5
		33.75
		50
		z, xi

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

	Задача 4_2.4. На отрезке [-l,l] задана функция y = f (x) .			Выполнить приближение этой функции при
		a0	m
	помощи тригонометрических многочленов Tm (x) =		+ åak cos(π kx) + bk sin(		π kx) .
	2		k=1	l	l

ПОРЯДОК РЕШЕНИЯ ЗАДАЧИ:

1.Задать степень тригонометрического многочлена m (выполнить задание для пяти значений m = 1,3,5,10 и 20.

2.Выполнить преобразование Фурье. Найти коэффициенты Фурье до m включительно.

ì			1		l
ïak	=
			l −òl
ï
í		1 l
ï
ïbk	=				ò
			l
î					−l

f (u) cos(		kπ		u)du,	k = 0,1,2,3,...,m
		l
f (u)sin(	kπ		u)du,		k = 1,2,3,...,m .
		l

3.Составить тригонометрический многочлен степени m

		a0	m
	Tm (x) =		+ åak cos(π kx) + bk sin(		π kx) .
	2		k=1	l	l

4.Вывести на один чертеж график функции y = f (x) на отрезке [-l,l] и графики всех тригонометрических многочленов, полученных для различных значений m. Сделать выводы.

Исходные данные:

Задана функция на промежутке:

f(x) := x + e	(− x)2			x := -2, -2 + 0.01.. 2					L := 2	f(-2) = 52.5981500331442
Решение задачи:										f(2) = 56.5981500331442
Koeff_a(m,L, f) :=				for	kÎ 0.. m
					1		óL	æ p	ö
				ak			ô
					¬ L	×	ô	f(x) × cosè L	× k × xø dx
							õ− L
				a
Koeff_b(m,L, f) :=				for	kÎ 0.. m
					1		óL	æ p	ö
				bk			ô
					¬ L	×	ô	f(x) × sinè L	× k × xø dx
							õ− L
				b

Составим функцию вычисляющую тригонометрический многочлен:

	a0		m	æ	æ p	ö	æ p	öö
T(x, m, L,a, b) :=			+ å	èak	× cosè L	× k × xø	+ bk × sinè L	× k × xøø
	2
			k = 1

Найдем тригонометрические многочлены степени m=1,3,5,10,20 и построим их графики:

m := 1

a := Koeff_a(m,L, f)			b := Koeff_b(m, L,f)
F_1(x) := T(x,m,L, a, b)
								100
	f(x)							50
	F_1(x)
							2		50		0			2
										x

m := 3
a := Koeff_a(m,L, f)	b := Koeff_b(m, L,f)
F_3(x) := T(x,m,L, a, b)

f(x)

F_3(x)

m := 5

a := Koeff_a(m,L, f) F_5(x) := T(x,m,L, a, b)

100

50

2

0

2

50

x

b := Koeff_b(m, L,f)

	100
f(x)	50
F_5(x)
2	0	2
	50
	x
m := 10
a := Koeff_a(m,L, f) b := Koeff_b(m, L,f)
F_10(x) := T(x,m,L, a,b)
	100
f(x)
F_10(x)	50
2	0	2
	x

m := 20

a := Koeff_a(m,L, f)

b := Koeff_b(m, L,f)

F_20(x) := T(x,m,L, a,b)

100

f(x)

	F_20(x)		50
		2	0		2
			x
Выведем на один график функцию y=f(x) на отрезке -2 до 2 и графики всех
тригонометрических многочленов, полученных при различных значениях m=1,3,5,10,20
f(x)					40
F_1(x)
F_3(x)
F_5(x)
F_10(x)					20
F_20(x)
	2	1.5	1	0.5	0	0.5	1	1.5	2
					x
Вывод:	При наблюдении за тригонометрическими многочленами, можно убедиться в
	том, что чем больше степень многочлена, тем точнее он приближается к
	исходной функции.

Задача 4_2.5. Функция y=f(t) , представляющая некоторый периодический сигнал, задана таблицей

значений y	i	= f (x )	на отрезке [0,2π ] в точках x	i	=	2π	×i , где N задает частоту дискретизации,
		i				N
i = 0,..., N -1. Известно, что в каждом из представленных значений уровней сигнала присутствует

некоторый шум, имеющий случайное распределение, т.е. на самом деле yi = f (xi ) + Noisei . Параметры

шума заранее не известны.

При помощи использования дискретного преобразования Фурье провести фильтрацию (аппроксимацию)

	сигнала,	включив в аппроксимирующий (сглаживающий) тригонометрический многочлен
		A0	m
	Qm (t) =		+ å Ak cos(kt) + Bk sin(kt) только информативные гармоники (по значениям амплитуд)
		2
			k =1

для соответствующей степени m. Вывести на один чертеж точечный график исходного «зашумленного» дискретного сигнала и непрерывный график сглаженного сигнала Qm (t) .

Теоретический материал:

Дискретным преобразованием Фурье называется действие, с помощью которого по заданной в интервале функции f (x) строится система чисел {ak ,bk }∞k =0 - называемых коэффициентами Фурье. По традиции,

именно эти числа также обозначают словами «преобразование Фурье» данной функции. При этом числа ak

называются косинус-преобразованием Фурье функции f (x) , а числа bk называются синус-

преобразованием Фурье функции f (x) . У преобразования Фурье функции есть множество свойств и применений; в частности, известен ряд прикладных вопросов, в которых ту или иную информацию о функции f (x) удается получить только через ее преобразование Фурье, которое в таких ситуациях удается получить некоторыми косвенными средствами.

Рассмотрим следующий частный случай. Функция f (x) рассматривается на интервале [0.2π ] и

притом только в его отдельных точках

2π

xl = n l, l = 0,1,...,n - 1

при некотором заранее заданном и фиксированном числе n . Значения функции f (x) в этих точках считаются известными; обозначим

fl = f (xl ), l = 0,1,...,n - 1.

В равенстве для ряда Фурье

					a0	∞
				f (x) =		+ å(ak cos(kx) + bk sin(kx))
	положим x = xl . Получим			2		k=1
				∞
			a0
	fl	=		+ å(ak		cos(kxl ) + bk sin(kxl ))
	(1)	2		k =1

Проанализируем соотношение (1). Если произвольное целое неотрицательное число k разделить с остатком на число n , то получится соотношение k = pn + q , где для целых p,q имеются лишь следующие возможности:

p= 0,1,2,..., q = 0,1,2,...,n -1.

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

			A0	n−1
	(2)	fl =		+ å( Aq cos(qxl ) + Bq sin(qxl )) ,
			2
				q=1

где

∞	∞
A0 = a0 + 2åapn , Aq = åapn+q ,
p=1	p=0

∞

B0 = 0, Bq = åbpn+q , p=0

q = 1,2,...,n −1.

Отметим, что теперь все суммы в (2) - конечны. Известен следующий факт о тригонометрических

суммах:

для всех чисел q =1,2,3,...n −1 имеют место равенства

n−1

åcos(qxl ) = 0,

l =0

n−1

åsin(qxl ) = 0.

l=0

Если обе части соотношения (2) умножить на cos(rxl ) и затем просуммировать по l , то, с учетом только что сказанного, легко получить, что

			2	n−1
	(3а)	Ar =		å fl cos(rxl ) ;
			n l=0

а если обе части (2) умножить на sin(rxl ) и, с учетом того же утверждения о суммировании косинусов и синусов, получим соотношение

		2	n−1
(3б)	Br =		å fl sin(rxl ),
		n l=0
причем в (3а) и (3б)	r = 0,1,2,...,n −1. Числа Ar ,Br ,r = 0,1,2,...,n −1, называются			дискретным
преобразованием Фурье функции f (x) . Если в равенстве (2) заменить xl на произвольный				x , то оно из

точного станет приближенным. Его правую часть в этом случае называю тригонометрической интерполяцией функции f (x) .

Исходные данные:

d := (8.89168197002832 7.23030753944281 9.04159870460986 9.04604086793189 8.09279400592598 5.3735

N := 118
i := 0.. N - 1	y := dT

x :=	2 × p	× i Узлы дискретизации
i	N

Узлы дескретизации и значения в них:

					0																																						0
		0			0																										0			8.89168197002832
		1	0.053247333111691																												1			7.23030753944281
		2	0.106494666223383																												2			9.04159870460986
		3	0.159741999335074																												3			9.04604086793189
		4	0.212989332446766																												4			8.09279400592598
		5	0.266236665558457																												5			5.37354256495133
		6	0.319483998670148																												6			4.00280118161518
x =		7	0.37273133178184																								y =				7			2.16299559425184
		8	0.425978664893531																												8			2.7334742563888
		9	0.479225998005223																												9			2.65644591191931
		10	0.532473331116914																												10			3.04700220933884
		11	0.585720664228606																												11			4.05828797365486
		12	0.638967997340297																												12			6.14701321256455
		13	0.692215330451988																												13			5.37817286193427
		14	0.74546266356368																												14			6.04807600866672
		15	0.798709996675371																												15			6.80259522317658
		yi	10
			0
			0		1															2										3				4									5																			6								7
																																				xi

1. Рассчитывам коэффициенты Фурье на основании известных формул Aq и Bq, q=0,1..,N-1:

A := for q Î 0.. N - 1

B := for q Î 0.. N - 1

		2	éN−1	ù		2	éN−1	ù
	Aq ¬		× êêå	(yl × cos(q × xl))úú	Bq ¬		× êêå	(yl × sin(q × xl))úú
		N				N
			ël = 0	û			ël = 0	û
	A				B

2. Находим амплитуды соответствующих гармоник:

q := 1.. N - 1 Amq := (Aq)2 + (Bq)2

	Amq		4
			2
			0
				0	6			12		18					24			30			36		42					48	54			60	66		72		78	84			90		96				102				108				114				120
																																q

3. Анализ гармоник (достаточно анализировать гармоники с номерами от 1 до N/2 - остальные симметричны первым) показывает, что основной полезный сигнал содержится в 8-ой и 10-ой гармониках (их амплитуды значительно превосходят амплитуды остальных гармоник). На основании этой информации выбираем оптимальную степень тригонометрического многочлена приближения (m=10).

m := 10

4. Соcтавляем тригонометрический многочлен степени m. Гармоники, имеющие порядок выше, чем m, мы отбрасываем, осуществляя тем самым высокочастотную фильтрацию.

			A20 +	m	(Aq × cos(q × t) + Bq × sin(q × t))
t := 0, 0 + 0.001.. 2p		Q(t) :=		å
				q = 1
	10
Q(t)
yi	0
	0	1	2		3	4	5	6	7
						t , xi

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

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

4256495133 4.00280118161518 2.16299559425184 2.7334742563888 2.65644591191931 3.04700220933884 4.