Математические методы проектирования

«Математические методы проектирования»

1. Моделирование сигнала на входе системы

2. Моделирование работы системы при входном сигнале

3. Моделирование шума с заданной плотностью распределения вероятностей

4. Моделирование работы системы при входном шуме Выводы Приложения

В данной курсовой работе требуется:

— выполнить моделирование входного заданного сигнала, построить график, амплитудный и фазовый спектр сигнала;

— построить математическую модель системы (рис. 2.1), рассчитать и построить амплитудно-частотные характеристики (АЧХ) каждого из блоков и всего устройства в целом, а также выполнить моделирование работы системы при заданном входном сигнале;

— выполнить моделирование шума с законом распределения вероятностей Рэлея и дисперсией D=12, оценить среднее значение и дисперсию отсчетов шума, построить гистограмму и проверить адекватность модели шума по критерию Пирсона;

— выполнить моделирование работы системы при входном шуме, построить графики спектра, автокорреляционной функции и гистограммы выходного сигнала, оценить статистические параметры и выбрать статистическую модель для выходного шума.

1. МОДЕЛИРОВАНИЕ СИГНАЛА НА ВХОДЕ СИСТЕМЫ

На входе моделируемой системы имеется гармонический сигнал:

s (t)=(1+аcos2рЩt)•sin (2рf0t)+sin2рaЩt. (1.1)

Необходимо выполнить тригонометрические преобразования и привести выражение (1.1) к виду, чтобы наглядно можно было определить теоретическое количество спектральных составляющих, их частоты и амплитуды. С учетом заданных значений a=0.8, Щ=50 Гц, f0=220 Гц получаем:

S (t)=sin (2*р*220*t)+0.4sin (2*р*170)+0.4sin (2*р*t*270)+sin (2*р*0.4*50*t).(1.2)

После преобразований стало видно, что сигнал состоит из четырех гармоник с амплитудами 1, 0.4, 1, 0.4 соответственно на частотах 50 Гц, 170 Гц, 220 Гц, 270 Гц. Для моделирования дискретного сигнала применялся следующий алгоритм:

— определяется максимальная частота, в данном случае fв=270 Гц;

— применяя теорему Котельникова, была рассчитана частота дискретизации

fd=2*fв 5=2700 Гц,(1.3)

где 5 — коэффициент запаса;

— определяем интервал времени между отсчетами:

Дt=1/fd =3,7· 10−4 с;(1.4)

— выбираем необходимое время моделирования:

T=1/fmin=0,1 с ,(1.5)

где fmin — минимальное значение частоты сигнала, в данном случае 10 Гц;

— определяем первоначальное значение отсчетов во временной области

N=T/Дt. (1.6)

Рассчитанный объем выборки равен N = 270 отсчетов. Однако, для использования в дальнейшем процедуры быстрого преобразования Фурье (БПФ) необходимо полученное значение N округлить до ближайшего верхнего числа, кратного N= 2k, где k — целое число. Соответственно принято значение N =512. В связи с этим требуется пересчитать шаг по времени и частоту дискретизации. Для N =512 шаг по времени Дt=1,9· 10−4 с.

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

Рисунок 1.1 -Моделируемый входной сигнал Подпрограмма для преобразований Фурье требует, чтобы поступающие в нее значения имели комплексный вид. Поэтому для получения комплексных значений входного сигнала была использована функция cmplx. Затем при помощи процедуры быстрого прямого преобразования Фурье (БППФ) был получен спектр сигнала. При использовании БППФ спектр урезается вдвое, а на интервале от N/2 до N спектр является комплексно сопряженным «зеркальным» отображением основной части. Как видно из уравнения входного сигнала, он имеет лишь мнимые составляющие. Действительная часть спектра равна 0.

Мнимая составляющая спектра (см. рис.1.2) вычисляется при помощи функции aimag.

Рисунок 1.2 — Мнимая составляющая спектра входного сигнала Амплитудный спектр (рис. 1.3) получают взятием модуля от комплексного сигнала. В нашем случае использовалась функция cabs.

Рисунок 1.3 — Амплитудный спектр входного сигнала Отношение мнимой составляющей спектра к действительной дает информацию о фазе данного процесса. В среде Fortran, была использована функция atan2. График фазового спектра входного сигнала изображен на рис. 1.4.

Рисунок 1.4 — Фазовый спектр входного сигнала

2. МОДЕЛИРОВАНИЕ РАБОТЫ СИСТЕМЫ ПРИ ВХОДНОМ СИГНАЛЕ

Цифровое моделирование работы системы заключается в последовательном моделировании преобразования сигнала по структурной схеме при известных передаточных характеристиках отдельных блоков. Структурная схема системы представлена четырьмя блоками (см. рис.2.1)

Рисунок 2.1 — Моделируемая система Схемы электрические принципиальные для блоков 9, 12, 14 приведены на рис. 2.2.

а) б) в) а) схема электрическая принципиальная блока 9;

б) схема электрическая принципиальная блока 12;

в) схема электрическая принципиальная блока 14.

Рисунок 2.2 -Схемы электрические принципиальные блоков

моделируемой системы Передаточные функции блоков 9, 12, 14 определяются следующим образом:

(2.1)

где p — оператор Лапласа;

Li — индуктивность;

Ri — сопротивление;

Сi — емкость;

Графики АЧХ блоков 9, 12, 14 изображены на рисунках 2.3 — 2.5.

Рисунок 2.3 — Амплитудно-частотная характеристика блока 9

Рисунок 2.4 — Амплитудно-частотная характеристика блока 12

Рисунок 2.5 — Амплитудно-частотная характеристика блока 14

Характеристики блока Т заданы таблично (см. таблицу 2.1).

Таблица 2.1

Для получения промежуточных точек таблицы использовался полином Лагранжа:

(2.4)

где

N — степень полинома Лагранжа;

fi — значения функции в узлах интерполяции;

x*- точка, в которой вычисляется значение полинома;

xi — значения в узлах интерполяции.

Графики амплитудно-частотной и фазо-частотной характеристик блока Т, полученных с помощью интерполяционной формулы Лагранжа, отображены на рис. 2.6 — 2.7.

Рисунок 2.6 — Амплитудно-частотная характеристика блока Т Рисунок 2.7 — Фазо-частотная характеристика блока Т Общая передаточная функция рассчитывается последовательным блоком операций с учетом правил для соединений блоков. В данной системе присутствуют элементы, соединенные с помощью отрицательной обратной связи (Wт и W14). После замены блоков, соединенных отрицательной обратной связью на единый блок, получаем последовательное соединение элементов системы. Общая передаточная функция такой системы равна произведению всех ее элементов. В результате общая передаточная функция примет вид:

(2.5)

где Wi — передаточная функция i-го блока.

График АЧХ системы изображен на рисунке 2.8.

Рисунок 2.8 — Амплитудно-частотная характеристика системы Алгоритм моделирования выходного сигнала состоит в применении быстрого обратного преобразования Фурье (БОПФ) к выражению:

(2.6)

где W (p) — передаточная функция системы в комплексной форме;

X (p) — спектр входного сигнала.

В результате моделирования был получен амплитудный спектр (рис. 2.9) и фазовый спектр (рис. 2.10) сигнала на выходе системы.

Рисунок 2.9 — Амплитудный спектр выходного сигнала Рисунок 2.10 — Фазовый спектр выходного сигнала Применяя обратное преобразование Фурье к уравнению (2.6) был получен сигнал на выходе системы. Для сравнения были построены входной и выходной сигнал (рис. 2.11). Так как амплитуда выходного сигнала значительно меньше амплитуды входного сигнала, то для графика, значения амплитуды выходного сигнала брались с коэффициентом 100.

Рисунок 2.11 — Входной и выходной сигналы При сравнении графиков входного и выходного сигналов (рис. 2.11), а также спектров входного сигнала (1.3 — 1.4) со спектрами выходного сигнала (2.9 — 2.10) можно заметить значительное влияние моделируемой системы на входной сигнал и его спектр.

3. МОДЕЛИРОВАНИЕ ШУМА С ЗАДАННОЙ ПЛОТНОСТЬЮ РАСПРЕДЕЛЕНИЯ ВЕРОЯТНОСТЕЙ

Требуется выполнить моделирование шума с законом распределения вероятностей Рэлея и дисперсией D=12, где у=.

Для получения реализаций шума с заданным законом распределения используется метод обратной функции. Чтобы получить алгоритм моделирования, приравниваются функция распределения и равномерная случайная величина, где r — сформированная при помощи мультипликативного алгоритма равномерная случайная величина [0.1]. После вычислений и преобразований был получен алгоритм моделирования вида, который позволил сформировать массив значений шума.

График моделируемого шума изображен на рисунке 3.1.

Рисунок 3.2 — Моделируемый шум Выполняя статистический анализ полученного случайного процесса, была определена оценка математического ожидания:

(3.4)

где xi — значения отсчетов шума;

N — количество отсчетов.

Оценка дисперсии:

.(3.5)

моделирование сигнал шум амплитудный

По результатам расчетов были получены следующие значения: =0,23 584,= 0,149 634.

Далее с помощью подпрограммы Rhist была построена гистограмма, показывающая количество попаданий случайных чисел в определённый интервал. Было взято 11 интервалов. Гистограмма представлена на рис. 3.3.

Рисунок 3.3 — Гистограмма входного шума В качестве проверки соответствия теоретического и экспериментального распределений используется критерий Пирсона. Требуется определить ч2и оценить степень соответствия распределений.

где M — количество интервалов гистограммы;

N — объем выборки;

Pk — площадь k-го интервала гистограммы, определяемая по методу Симпсона;

nk — экспериментально полученное количество попаданий в k-й интервал гистограммы.

Был построен график сравнения теоретического и экспериментального количества попаданий в определенный интервал (рис 3.4).

Рисунок 3.4 — Сравнение теоретического и экспериментального числа попаданий в определенный интервал По итогам расчета имеем ч2 = 4,922. Вероятность соответствия теоретического закона и эмпирического для полученного ч2 является 0,8.

4. МОДЕЛИРОВАНИЕ РАБОТЫ СИСТЕМЫ ПРИ ВХОДНОМ ШУМЕ

Используя формулу (2.6), может быть вычислен спектр выходного шума. В результате был получен амплитудный спектр шума на выходе системы (см. рис.4.1).

Рисунок 4.1 — Амплитудный спектр выходного шума Можно заметить, что спектр выходного шума по форме похож на передаточную характеристику заданной системы.

Выходной сигнал был восстановлен из спектра при помощи обратного преобразования Фурье (см. рис. 4.2).

Рисунок 4.2 — Шум на выходе системы С помощью теоремы Винера-Хинчина где — обратное преобразование Фурье;

S2(f) — спектральная плотность мощности, была получена автокорреляционная функция (АКФ) (см. рис. 4.3).

Рисунок 4.3 — Автокорреляционная функция Гистограмму для выходного шума строим по алгоритму, по которому строилась гистограмма выходного сигнала (см. рис. 4.4).

Рисунок 4.4 — Гистограмма выходного шума Чтобы узнать, по какому закону распределен выходной сигнал, необходимо определить эксцесс и асимметрию .

(4.2)

где ni — отсчет выходного шума;

— оценка математического ожидания;

— оценка среднеквадратического отклонения Эксцесс определяется как:

.(4.3)

Получили такие результаты: =7,946· 10−5, = 3,1 942 984· 10−4, =0,0779, = - 0,47. Наиболее близким к полученному распределению оказался нормальный закон, у которого асимметрия и эксцесс равны 0 (см. рис. 4.5)

После проверки сходства законов по критерию Пирсона, получаем ч2=28,688. Вероятность соответствия теоретического закона и эмпирического для полученного ч2 является <0,05

График сравнения теоретического и экспериментального числа попаданий в определенный интервал представлен на рис 4.6.

Рисунок 4.5 — Фрагмент плоскости моментов Рисунок 4.6 — Сравнение теоретического и экспериментального числа попаданий в определенный интервал

ВЫВОДЫ

В работе было выполнялось моделирование гармонического сигнала. По результатам расчетов спектр сигнала на входе имеет четыре гармонических составляющих. Глядя на полученный график (рис. 1.3) можно отметить, что экспериментально полученный амплитудный спектр совпадает с теоретически ожидаемым.

Что касается передаточной функции, из формулы (2.6) следует, что на сигнал и его спектр значительным образом влияют параметры системы. Сравнение графиков спектров сигнала на входе устройства (см. рис.1.3 — 1.4) и на выходе (см. рис.2.9 — 2.10) показывает, что наблюдается искажение спектра.

Было выполнено моделирование шума с распределением вероятностей Рэлея. Совмещенный график теоретического распределения и гистограммы иллюстрирует достаточно высокую степень согласия. По итогам расчета имеем ??2 = 4,92. Экспериментальные данные соответствуют теоретическому распределению с вероятностью 0,8.

При прохождении через систему, ее воздействие существенно искажает форму шума. Поэтому для определения закона распределения выходного шума необходимо рассчитать эксцесс и асимметрию. Полученные значения были нанесены на плоскость моментов (рис. 2.5). С их помощью был определен ближайший закон распределения. Наиболее оптимальным оказался нормальный закон распределения. После проверки соответствия экспериментально полученного закона распределения и нормального распределения по критерию Пирсона оказалось, что вероятность соответствия менее 0,05.

Приложение А

Моделирование заданного входного сигнала

Program spectr

integer, parameter: N=512

real a/0.8/, Om/50/, f0/220/

real T/0.1/, pi/3.1 415 927/, fmin/10/

real t_mas (N), S (N), Spectr_Re (N), Spectr_Im (N), Sp_A (N), Sp_Ph (N), F (N)

integer i

complex cS (N)

dt=T/real (N)

do i=1,N

t_mas (i)=real (i-1)*dt

S (i)=sin (2*pi*220*t_mas (i))+0.4*(sin (2*pi*t_mas (i)*170)+sin (2*pi*t_mas (i)*270))+sin (2*pi*0.4*50*t_mas (i))

end do

open (unit=1,file='Время.txt', action='WRITE')

write (1,*) t_mas

close (1)

do i=1,N/2

F (i)=real (i-1)*fmin

end do

call Draw (N, 1, t_mas, S, Red, 'Line', 'signal')

accept *

do i=1,N

cS (i)=cmplx (S (i), 0.0)

end do

open (unit=1,file='vhod.txt', action='WRITE')

write (1,*) cS

close (1)

call FFT (cS, N,0)

do i=1,N/2

Spectr_Re (i)=real (cS (i))

Spectr_Im (i)=aimag (cS (i))

end do

call Draw (N/2, 1, F, Spectr_Re, 4, 'Stick', 're spectr')

accept *

call Draw (N/2, 1, F, Spectr_Im, 4, 'Stick', 'im spectr')

accept *

do i=1,N/2

Sp_A (i)=cabs (cS (i))

if ((abs (real (cS (i)))<1.0e-5).and.(abs (aimag (cS (i)))<1.0e-5)) then

Sp_Ph (i)=0

else

Sp_Ph (i)=atan2(aimag (cS (i)), real (cS (i)))

end if

end do

call Draw (N/2, 1, F, Sp_A, 0, 'Stick', 'амплитудный спектр')

accept *

call Draw (N/2, 1, F, Sp_Ph, 0, 'Stick', 'фазовый спектр')

end program spectr

Приложение Б

Моделирование работы системы при входном сигнале Расчет характеристик блока Т

program first

implicit none

parameter N=11,m=301

real f_z (N)/0,100,200,300,400,500,600,700,800,900,1000/

real x_z (N)/-0.8,-2.05,-4.45,-6.78,-8.77,-10.4,-11.9,-13.1,-14.2,-15.2,-16.1/

real y_z (N)/0,-29.7,-48.8,-59.7,-66.4,-70.7,-73.7,-76,-77.7,-79,-80.0/

real ach (m), fch (m), polinom, f (m), t

integer i

t=(f_z (n)-f_z (1.0))/(m-1.0)

do i=1,m

f (i)=f_z (1.0)+t*(i-1.0)

end do

do i=1,m

ach (i)=polinom (f_z, x_z, N, f (i))

end do

open (unit=1,file='АЧХ.txt', action='WRITE')

write (1,*) ach

close (1)

do i=1,m

fch (i)=polinom (f_z, y_z, N, f (i))

end do

open (unit=2,file='ФЧХ.txt', action='WRITE')

write (2,*) fch

close (1)

call Draw (m, 1, f, ach, 3, 'Line', 'график АЧХ')

accept *

call Draw (m, 1, f, fch, 0, 'Line', 'график ФЧХ')

accept *

end program first

Расчет общей передаточной функции и моделирование выходного сигнала

implicit none

integer i

parameter N=512, R1=5.1, R2=0.1, C1=68.0e-6, L1=1.0e-2, M=N/2, Tn=0.1

real pi/3.1 415 926/, fmin/10.0/, dt, Sp_F (M), Sp_A (M)

real W9(N/2), W12(N/2), W14(N/2), Wt (N/2), W (N/2), f (N/2), ach (n/2), fch (N/2), Y (N), t (N, 2), s (N, 2)

complex cj/(0.0,1.0)/, p (N/2), cW9(N/2), cW12(N/2), cW14(N/2), cW (N), cWt (N/2), cs (N), cv (N)

do i=1,N/2

f (i)=real (i-1)*fmin

end do

open (unit=1,file='Частота.txt', action='WRITE')

write (1,*) f

close (1)

do i=1,N/2

p (i)=cj*2*pi*f (i)

cW9(i)=(p (i)*L1+R2)/(p (i)**2*R1*C1*L1+p (i)*(R1*R2*C1+L1)+R1+R2)

cW12(i)=R2/(p (i)*L1+R1+R2)

cW14(i)=(p (i)**2*C1*L1+p (i)*C1*R2)/(p (i)**2*C1*L1+p (i)*C1*(R1+R2)+1.0)

end do

W9=cabs (cW9)

W12=cabs (cW12)

W14=cabs (cW14)

open (unit=1,file='.АЧХ.txt', action='READ')

read (1,*) ach

close (1)

open (unit=2,file='.ФЧХ.txt', action='READ')

read (2,*) fch

close (2)

do i=1,n/2

ach (i)=10**(ach (i)/20)

fch (i)=pi*fch (i)/180

cWt (i)=ach (i)*exp (cj*fch (i))

end do

Wt=cabs (cWt)

call Draw (N/2, 1, f, W9, 0, 'Line', 'АЧХ 9')

accept *

call Draw (N/2, 1, f, W12, 1, 'Line', 'АЧХ 12')

accept *

call Draw (N/2, 1, f, W14, 2, 'Line', 'АЧХ 14')

accept *

call Draw (N/2, 1, f, Wt, 2, 'Line', 'АЧХ T')

accept *

cW (1:M)=cWt*cW9*cW12/(1+cWt*cW14)

W (1:M)=cabs (cW (1:M))

call Draw (N/2, 1, f, W, 4, 'Line', 'АЧХ общее')

accept *

do i=1, N/2−1

cW (N/2+i+1)=conjg (cW (N/2-i+1))

end do

open (unit=1,file='Передат.txt', action='WRITE')

write (1,*) cW

close (1)

open (unit=1,file='.secondvhod.txt', action='READ')

read (1,*) cS ! signal

close (1)

s (1:N, 1)=real (cs)

call FFT (cS, N,0) ! cs= input spectrum

cv=cS*cW ! output spectrum

do i=1,M

Sp_A (i)=abs (cv (i))

if ((abs (real (cv (i)))<1.0e-5).and.(abs (aimag (cv (i)))<1.0e-5)) then

Sp_F (i)=0

else

Sp_F (i)=atan2(aimag (cv (i)), real (cv (i)))

end if

end do

call FFT (cv, N,1) ! cv — output complex signal

call Draw (50, 1, f, Sp_A, 4, 'Stick', 'Амплитудный спектр')

accept *

call Draw (50, 1, f, Sp_F, 4, 'Stick', 'Фазовый спектр')

accept *

dt=Tn/real (N-1.0)

do i=1,N

t (i, 1)=real (i-1)*dt

t (i, 2)=t (i, 1)

Y (i)=real (cv (i)) ! output real signal

end do

open (unit=1,file='Время.txt', action='WRITE')

write (1,*) t

close (1)

s (1:N, 2)=100*Y

call Draw (N, 2, t, s, 4, 'Line', 'Выходной сигнал (t)')

end

Приложение В

Моделирование шума с законом распределения вероятностей Рэлея

program shumec

implicit none

parameter D=12, N=512, Int=11, M=21, Nf=int*30+1,nt=31

integer i, j, k

realsi, r/0.15 846 827/, Sh (N), Mx, s, Dx, SKO, Z (Int), Hist (Int), t_mas (N), Pl (Int), x (nf)

real rmin, rmax, simpson, ft (nf), Xi, X2(Int, 2), Y2(Int, 2)

si=sqrt (0.43/D)

do i=1,N

r=r*37

r=r-aint®

Sh (i)=sqrt (-2*si**2*log (1-r))

end do

open (unit=1,file='.secondВремя.txt', action='READ')

read (1,*) t_mas

close (1)

call Draw (N, 1, t_mas, Sh, 4, 'Line', 'Шум')

accept *

open (unit=1,file='Шум.txt', action='WRITE')

write (1,*) Sh

close (1)

s=0.0

do i=1,N

s=s+Sh (i)

end do

Mx=s/real (N)

print *,'Мат ожидание:', Mx

accept *

Dx=0.0

do i=1,n

Dx=Dx+((Sh (i)-Mx)**2)

end do

Dx=Dx/real (N-1)

Print *,'Дисперсия:', Dx

accept *

SKO=sqrt (Dx)

Print *,'СКО:', SKO

accept *

rmin=Sh (1);rmax=Sh (1)

do i=2,N

if (Sh (i)>rmax) rmax=Sh (i)

if (Sh (i)

end do

dx=(rmax-rmin)/(int-1)

do i=1,Int

Z (i)=rmin+dx*real (i-1)

end do

call Rhist (Hist, Int, Sh, N, rmin, rmax)

Call Draw (Int, 1, Z, Hist, 4, 'Box', 'Гистограмма')

accept *

dx=(rmax-rmin)/(nf-1)

do i=1,nf

x (i)=rmin+dx*(i-1)

ft (i)=(x (i)/si**2)*exp ((-x (i)**2)/(2*si**2))

end do

s=1

do i=1, Int

Pl (i)=simpson (ft (s), nt, x (s), x (s+30))

s=s+30

end do

s=sum (pl)

!print *,'x=', s

Xi=0

do i=1,Int

Xi=Xi+(Hist (i)-Pl (i)*N)**2/(Pl (i)*N)

end do

print *,'xi=', Xi

accept *

do i=1,Int

X2(i, 1)=z (i)

X2(i, 2)=z (i)

Y2(i, 1)=hist (i)

Y2(i, 2)=pl (i)*N

end do

call Draw (int, 2, X2, Y2, (/2,7/), (/'Box','Line'/), '1')

accept*

end program shumec

Приложение Г

Моделирование работы системы при входном сигнале шума.

program mmm

implicit none

parameter N=512, M=N/2, Int=11, Nf=int*30+1,nt=31

integer i

real Sh (N), f (M), s, Sh_v (M), AKF (N), Hist (Int), z (Int), Mx, SKO, alfa, eps, rmin, rmax, t (N), simpson

real Pl (int), Xi, dx, x (nf), ft (nf), X2(Int, 2), Y2(Int, 2)

complex cW (N), AKFc (N), Sp (N), cSh (N)

open (unit=1,file='.thirdПередат.txt', action='READ')

read (1,*) cW (1:N)

close (1)

open (unit=1,file='.thirdЧастота.txt', action='READ')

read (1,*) f

close (1)

open (unit=1,file='.fourthШум.txt', action='READ')

read (1,*) Sh

close (1)

cSh=cmplx (Sh)

call FFT (cSh, N,0)!spektr shuma

do i=1,N

cSh (i)=cSh (i)*cW (i)

end do

Sh_v (1:M)=abs (cSh (1:M))

call Draw (M, 1, f, Sh_v, 0, 'Stick', 'Спектр выходного шума')

accept *

open (unit=1,file='.thirdВремя.txt', action='READ')

read (1,*) t

close (1)

Call FFT (cSh, N, 1)!signal shum

Sh=real (cSh)

call Draw (N, 1, t, Sh, 0, 'Line','Shum')

accept *

mx=sum (Sh)/N

SKO=sqrt (sum ((Sh-mx)**2)/(n-1))

print *, mx, SKO

accept *

s=0.0

do i=1,N

s=s+(Sh (i)-Mx)**3

end do

alfa=s/(N*SKO**3)

print *,'Ассимеирия=', alfa

accept *

s=0.0

do i=1,N

s=s+(Sh (i)-Mx)**4

end do

eps=s/(N*SKO**4)-3.0

print *,'Эксцесс=', eps

accept *

rmin=Sh (1);rmax=Sh (1)

do i=2,N

if (Sh (i)>rmax) rmax=Sh (i)

if (Sh (i)

end do

do i=1,Int

Z (i)=rmin+((rmax-rmin)/(Int-1.0))*real (i-1.0)

end do

call Rhist (Hist, Int, Sh, N, rmin, rmax)

call Draw (Int, 1, Z, Hist, 0, 'Box', 'Гистограмма')

accept *

rmax=maxval (Sh)

rmin=minval (Sh)

dx=(rmax-rmin)/(nf-1.0)

do i=1,Int

Z (i)=rmin+dx*real (i-1.0)

end do

do i=1,nf

x (i)=rmin+dx*(i-1.0)

ft (i)=((1.0/(sqrt (2*3.1415*sko)))*exp (-(x (i)-mx)**2/(2*sko**2)))*sqrt (1/sko)

print *,'x=', x (i)

end do

call Draw (nf, 1, x, ft, 0, 'Line', 'zn')

accept *

s=1

do i=1, Int

Pl (i)=simpson (ft (s), nt, x (s), x (s+30))

s=s+30

end do

s=sum (pl)

print *,'x=', s

Xi=0

do i=1,Int

Xi=Xi+(Hist (i)-Pl (i)*N)**2/(Pl (i)*N)

end do

print *,'x=', Xi

accept *

do i=1,Int

X2(i, 1)=z (i)

X2(i, 2)=z (i)

Y2(i, 1)=hist (i)

Y2(i, 2)=pl (i)*N

end do

call Draw (int, 2, X2, Y2, (/2,7/), (/'Box','Line'/), '1')

accept*

do i=1,M

Sh (i)=Sh (i)-Mx

end do

cSh=cmplx (Sh)

call FFT (cSh, N,0)

do i=1,N

AKFc (i)=cSh (i)*conjg (cSh (i))

end do

call FFT (AKFc, N,1)

AKF (1:M)=real (AKFc (1:M))

AKF=AKF/AKF (1)

call Draw (N/2, 1, t, AKF, 0, 'Line', 'АКФ шума')

accept *

end program mmm

Приложение Д

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

Полином Лагранжа

real function polinom (x_z, f_z, N, x)

implicit none

integer i, k, N

real x_z (N), f_z (N), x, chisl, znam, w (N), sum

do k=1,N

chisl=1.0

znam=1.0

do i=1,n

if (i/=k)chisl=chisl*(x-x_z (i))

if (i/=k) znam=znam*(x_z (k)-x_z (i))

end do

w (k)=chisl/znam

end do

sum=0.0

do i=1,n

sum=sum+(w (i)*f_z (i))

end do

polinom=sum

end function polinom

Вычисление площади по формуле Симпсона

real function simpson (f, N, a, b)

integer N

real f (N), a, b, S_ch, S_nch, Pl

S_ch=0.0

do i=2,N-1,2

S_ch=S_ch+f (i)

end do

S_nch=0.0

do i=3,N-2,2

S_nch=S_nch+f (i)

end do

S=0.0

S=S+4*S_ch+2*S_nch+f (1)+f (N)

Pl=S*((b-a)/(3*(n-1)))

simpson=Pl

end function simpson