- •1. Теоретическое занятие
- •Компьютерная алгебра.
- •Использование систем компьютерной алгебры.
- •2. Практическое занятие
- •Лабораторная работа 1.
- •Разложение в ряд Тейлора
- •Лабораторная работа 2
- •Сходимость ряда Тейлора
- •Лабораторная работа 3
- •Вычисление интегралов
- •3. Практическое занятие
- •Лабораторная работа 4
- •Метод наименьших квадратов
- •Лабораторная работа 5
- •Быстрое преобразование Фурье
- •4 Практическое занятие
- •Задачи:
- •Литература
Таким образом, после измерения мы получим "экспериментальные" данные
>DataN:=[ [Rn[k],Jn[k]] $k=1..N]: A:=plot(DataN,style=point,symbol = circle,color=maroon): %;
По этим “экспериментальным”данным нам необходимо найти напряжение U. Для обработки данных по методу наименьших квадратов загрузим библиотеку
> with(CurveFitting):
и воспользуемся методом наименьших квадратов, используя команду LeastSquares, в которой мы явно указываем предполагаемый нами вид зависимости J=U/R
> g:=LeastSquares(DataN, r, curve=a/r); g := 103.0584229
r
Итак, что искомое напряжение примерно равно U=103. Проверьте, что увеличение числа измерений N приводит к улучшению результата.
Сравним график теоретической зависимости и экспериментальных данных
>B:=plot(g(r),r=Rn[1]..Rn[N],color=red): display([A,B]);
Заметим, что по этим данным мы можем построить не только гиперболу, но и произвольную функцию, например
>h:=LeastSquares(DataN, r, curve=a*exp(-0.6*r)+c);
>B:=plot(h(r),r=Rn[1]..Rn[N],color=red): display([A,B]);
Чтобы оценить, какую функцию надо использовать для аппроксимации по методу наименьших квадратов, мы воспользуемся процедурой
> describe[linearcorrelation]([Rn[k] $k=1..N],[Jn[j] $j=1..N]): evalf(%);
-0.7943085620
которая описывает линейные корреляции между данными (т.е. полиномиальные соответствия). В нашем случае мы получили, что наши данные лучше всего связываются законом
J =UR −0.79
Так как погрешность мы вводим генератором случайных чисел, то вы можете получить другой ответ.
Лабораторная работа 5
Быстрое преобразование Фурье
Рассмотрим функцию f(t) от одной переменной или сигнал, который зависит от времени.
Загрузим необходимые нам библиотеки процедур
11
> restart: with(linalg): with(plots):
и определим интервал [0,T] на котором мы собираемся изучать поведение этой функции
> T:=0.16;
Определим стационарный сигнал - сумму двух периодических функций, частоты которых
> w1:=50: w2:=200:
не меняются со временем
> f:=t-> sin(2*Pi*w1*t)+sin(2*Pi*w2*t); p1:=plot(f(t), t=0..T, title=`Сигнал 1`): %;
Используем метод дискретизации, при котором сигнал f(t) заменяется совокупностью его мгновенных значений, называемых выборками, или
отсчетами (samples).
Шаг дискретизации должен быть таким, чтобы было возможно восстановление непрерывной функции по ее отсчетам с допустимой точностью. По теореме Шенона-Котельникова этот шаг не должен превышать значения
> dT:=evalf(Pi/max(w1,w2));
Здесь max(w1,w2) - максимальное значение частоты в спектре сигнала.
Далее мы добавим к этоиу сигналу высокочастотный "белый" шум и, поэтому, заранее выберем шаг дискретизации меньше порогового значения
> dT:=0.001;
Зададим также величину выборки - т.е. число отсчетов N
> m:=10: N:=2^m;
Теперь все готово для того, чтобы перейти от непрерывной функции f(t) к дискретному набору ее значений. Определим массивы данных
>fk:= array(1..N): x:= array(1..N):
и заполним их
>for i to N do
fk[i]:=evalf(f(i*dT)); x[i]:=i*dT; end do:
Добавим теперь к этим значениям произвольный шум
> r:=randvector(N): fk:=evalm(fk+.02*r):
и сравним полученный нами "реальный" сигнал с "идеальным" сигналом
>maxN:=T/dT: Signal:=[ [x[j],fk[j]] $j=1..maxN]:
>pS:=plot(Signal,color=COLOR(RGB, 0.196, 0.6, 0.8), title="Сигнал + Шум"):
display([pS,p1]);
Естественно возникают два вопроса - какую информацию об исходной функции можно получить из "реального" сигнала и как очистить этот сигнал от шума?
Для решения этой задачи применяют быстрое преобразование Фурье
>im:=vector(N, 0):
>FFT(m,fk,im);
12
преобразование Фурье является комплексным преобразованием и, поэтому, нам необходимо задать мнимую составляющую сигнала - im, которая равна нулю в данном случае.
Вычислим максимальную частоту
> dw:=2*Pi/(N*dT);
отвечающую используемой нами выборке и определим ось частот
> fw := x -> (x-1)*evalf(dw/(2*Pi)): w := vector(N, fw):
для того, чтобы построить спектр сигнала
> Data:=[ [w[j],abs(fk[j])] $j=1..N]:
> p2:=plot(Data,color=coral, title="Спектр"):
%;
Ярко выраженные максимумы при ω1 = 50 и ω2 = 200 указывают на наличие в сигнале двух основных частот.
Используя эту информацию построим фильтр - маску, шириной
> dW:=9:
Фильтр - это функция
> filtr:=k->piecewise(k>w1-dW and k<w1+dW or k>N-w1-dW and k<N- w1+dW,1, k>w2-dW and k<w2+dW or k>N-w2-dW and k<N-w2+dW,1,0):
которая выглядит следующим образом
> plot(filtr,0..N,color=cyan, thickness=2);
Наложим этот фильтр на спектр сигнала
> for i from 1 to N do fk1[i]:=filtr(i)*fk[i]: im1[i]:=filtr(i)*im[i]: end do:
и посмотрим на полученный результат
> plot([[w[j],abs(fk1[j])] $j=1..N],color=coral, title="Спектр+Фильтр");
Теперь восстановим сигнал с помощью обратного преобразования Фурье
> iFFT(m,fk1,im1);
и сравним очищенный от шумов результат с исходной непрерывной функцией
> SF:=[ [x[j],fk1[j]] $j=1..maxN]: pF:=plot(SF,color=COLOR(RGB, 0.2, 0.5, 0.4), title="Очищенный сигнал"): display([pS,p1]);
display([pF,p1]);
Видно, что мы добились значительной очистки шумов. Степень очистки сигнала от шума можно регулировать, используя ширину маски. Попробуйте изменить ширину маски и посмотреть, как она влияет на степень очистки сигнала от шумов.
13