Методы цифрового моделирования

Размещено на http://www.allbest.ru/

ЗАДАНИЕ

На выполнение курсовой работы

«Математические методы и модели конструкторско-технологического проектирования»

1. Выполнить моделирование сигнала s(t) = (1 + a•cos2рЩt)•(sinрf0 + sin2рf0), построить s(t), спектр сигнала, где: а = 0.75, Щ = 20 Гц, f0 = 150 Гц.

2. Построить математическую модель системы (см. рис. 2.1), рассчитать и построить графики амплитудно-частотных характеристик каждого из блоков и всего устройства в целом. ЧХ блока «Т» задана таблично (см табл. 1). Выполнить моделирование работы системы при входном сигнале s(t). Построить графики выходного сигнала и его спектра.

Таблица 1.1 - АЧХ и ФЧХ блока «Т»

Параметры	Значения
F, Гц	0	50	100	150	200	250	300	350	400	450	500
АЧХ, дБ	0	-6,54	-11,5	-15,8	-18,3	-19,9	-21,7	-24	-25,2	-28,2	-32,1
ФЧХ, град	0	-58	-72,7	-78,3	-81,2	-83	-84,2	-85	-85,7	-86,2	-86,6

3. Выполнить моделирование шума n(t) с экспоненциальным законом распределения вероятностей с математическим ожиданием mn = 0.3. Оценить среднее значение и дисперсию отсчетов n(t), построить гистограмму и проверить адекватность модели n(t) по критерию Пирсона.

4. Выполнить моделирование работы системы при входном сигнале n(t). Построить графики спектра, автокорреляционной функции и гистограммы выходного сигнала. Выбрать статистическую модель для выходного сигнала, найти оценки ее параметров.

5. Оформить расчетно-пояснительную записку согласно ДСТУ3008-95.

математический цифровой моделирование

ВВЕДЕНИЕ

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

Объекты исследования: детерминированный сигнал; математическая модель системы; шум, распределенный по экспоненциальному закону; работа системы при воздействии на нее шума.

Метод исследования - статистическое моделирование в системе Fortran.

Ожидаемые результаты: моделирование всех объектов исследования; получение графиков сигналов и процессов; выбор статистической модели системы.

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

1.1 Расчет параметров моделирования

Задан сигнал s(t) = (1 + a • cos2рЩt) • (sinрf0 + sin2рf0). Раскрыв скобки и применив тригонометрическую формулу: sinA • cosB = sin(A - B) + sin(A + B), получаем

s(t) = sin2(рf0/2)t + sin2рf0t + sin2р()t + sin2р()t +

+ sin2р(f0 - Щ)t + sin2р(f0 + Щ)t. (1.1)

Подставим численные значения параметров а, Щ, f0 и получим

s(t) = sin2р75t + sin2р150t + 0,375sin2р55t + 0,375sin2р95t + 0,375sin2р130t +

+0,375sin2р170t(1.2)

Было выбрано наибольшее значение частоты - 170 Гц. Теперь используется теорема Котельникова, которая гласит, что, если аналоговый сигнал x(t) имеет ограниченный спектр, то он может быть восстановлен однозначно и без потерь по своим дискретным отсчётам, взятым с частотой более удвоенной максимальной частоты спектра Fmax:

,(1.3)

где Fmax -- верхняя частота в спектре, или (формулируя по-другому) по отсчётам, взятым с периодом чаще полупериода максимальной частоты спектра;

k - коэффициент запаса и k = 20.

Время дискретизации ?t рассчитывается по формуле:

?t=1/ Fd.(1.2)

После расчета получаем ?t = 2,9412 • 10-4 с.

Время наблюдения процесса рассчитывается так:

,(1.3)

где fмин - минимальный шаг по частоте, fмин = 5 Гц. Т = 0,2 с.

Количество отсчетов N (объем выборки) связано с необходимым временем наблюдения T за процессом:

N=T/?t.(1.4)

Для анализа спектрального состава процессов используют, как правило, процедуру быстрого преобразования Фурье (БПФ), которая в Fortran обозначается как FFT(*). Существует прямое и обратное преобразование Фурье. С помощью прямого получают спектр, а с помощью обратного сигнал восстанавливают из спектра. Особенностью процедуры БПФ является требование к объему выборки N= 2м, где М - целое число. Поэтому количество отсчетов N берем равное 1024. После этого необходимо пересчитать время дискретизации по времени ?t, которое равно 1,9531 • 10-4 с.

1.2 Описание алгоритма

Для выделения мнимой составляющей в Фортране используется функция aimag, для выделения реальной - функция real. Амплитудный спектр получают с использованием cabs - взятие модуля из комплексного числа. Фазовый спектр получают с заданной точностью 10-3, для его определения используется функция atan2, которая берется от комплексного массива сигнала. Для построения всех спектров используется М = N/2 точек, т. к. спектры обладают свойством зеркальности.

В результате моделирования получены следующие графики:

а) сигнал на входе системы;

б) амплитудный спектр входного сигнала;

в) фазовый спектр входного сигнала;

г) мнимая составляющая входного сигнала;

д) реальная составляющая входного сигнала.

Рисунок 1.1 - Моделируемый входной сигнал

Рисунок 2.2 - Амплитудно-частотная характеристика входного сигнала



Рисунок 2.3 - Фазо-частотная характеристика входного сигнала

Рисунок 2.4 - Мнимая составляющая спектра входного сигнала

Рисунок 2.5 - Реальная составляющая спектра входного сигнала

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

2.1 Описание математической модели системы и ее составляющих

На рисунке 2.1 изображена блок-схема заданной системы:

Рисунок 2.1 - Блок-схема системы

Составляющие систему блоки имеют следующий вид:

а) - схема блока 9

б) - схема блока 11

в) - схема блока 18

Рисунок 2.2 - Блоки 9, 11, 18

Математическое описание блоков - это их передаточные функции:

а) для блока 9

,(2.1)

б) для блока 11

, (2.2)

в) для блока 18

,(2.3)

где = j•щ - оператор Лапласа;

L - индуктивность;

R - сопротивление;

С - емкость.

Для блока Т АЧХ и ФЧХ заданы таблично (см. табл. 1.1). Для построения этих характеристик используется полином Лагранжа

,(2.4)

где хi - значения в узлах интерполяции.

Для получения графиков АЧХ и ФЧХ Т-блока используется подпрограмма интерполяции.

Общая передаточная функция была рассчитана в несколько этапов с применением правил для соединений блоков:

а) блок 11 и Т-блок соединены последовательно - их передаточные характеристики перемножаются;

W1 = W11•WT,(2.5)

б) новый полученный блок № 1 и блок № 18 - это цепь с отрицательной обратной связью;

W2 = (2.6)

в) новый полученный блок № 2 и блок № 9 соединены последовательно - их передаточные характеристики перемножаются.

Wобщ = W2•W9,(2.7)

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

2.2 Алгоритм моделирования

Для того, чтобы получить спектр выходного сигнала нужно

,(2.4)

где W(p) - передаточная функция;

X(p) - спектр воздействия, найден с помощью прямого преобразования Фурье.

Для того, чтобы получить выходной сигнал, используется обратное преобразование Фурье. Сигнал получается на М точек, а нужно на N, поэтому достраивается “зеркальная часть” передаточной характеристики с использованием функции conjg - комплексное сопряжение.

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

а) передаточные характеристики всех блоков (9, 11, 18) (см. рис. 2.6 - 2.8);

б) интерполированные АЧХ и ФЧХ Т-блока (см. рис. 2.9 - 2.10);

в) общая передаточная характеристика (см. рис. 2.11);

г) сигнал на выходе системы (см. рис. 2.12);

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

Рисунок 2.6 - Передаточная характеристика блока 9

Рисунок 2.7 - Передаточная характеристика блока 11



Рисунок 2.8 - Передаточная характеристика блока 18

Рисунок 2.9 - Интерполированная АЧХ Т-блока

Рисунок 2.9 - Интерполированная ФЧХ Т-блока

Рисунок 2.10 - Общая передаточная характеристика



Рисунок 2.11 - Сигнал на выходе системы

Рисунок 2.12 - Амплитудный спектр выходного сигнала

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

Рисунок 2.14 - Мнимая составляющая выходного сигнала



Рисунок 2.15 - Реальная составляющая выходного сигнала

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

Заданный шум распределен по экспоненциальному закону.

n(t) = л•e-л•t,(3.1)

где л - параметр.

Для моделирования задано математическое ожидание МХ = 0,3. Параметр л связан с математическим ожиданием

МХ = 1/ л.(3.2)

Из (3.2) был найден параметр л = 3,333.

Для имитации реализаций непрерывных случайных величин применяют метод обратной функции, позволяющий получить из реализаций ri - стандартной случайной величины реализации yi случайной величины с требуемым законом распределения. Для вычисления стандартной случайной величины был применен мультипликативный алгоритм.

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

Xi = - .(3.2)

Шум имеет вид, показанный на рисунке 3.1.



Рисунок 3.1 - Шум, распределенный по экспоненциальному закону

Для полученных непрерывных случайных величин рассчитана оценка дисперсии

,(3.3)

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

,(3.4)

Результаты расчета: = 9,2 • 10-2, у = 0,3.

Для построения гистограммы (см. рис. 3.2) используется подпрограмма Rhist.

Выбрано для построения 15 интервалов.



Рисунок 3.2 - Гистограмма входного шума

Для проверки согласия теоретического и экспериментального распределений используется критерий Пирсона. По значению ч2 и количеству степеней свободы ф = 13 определяется вероятность того, на сколько согласуются вышеуказанные распределения.

,(3.5)

где NPk - теоретическое количество попаданий в интервал гистограммы;

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

Для того, чтобы получить теоретическое количество попаданий в n-й интервал гистограммы, были посчитаны площади каждого интервала с использованием формулы Симпсона. Далее площадь каждого интервала была умножена на количество отсчетов (N = 1024).

После расчета получаем ч2 = 11,68 и вероятность соответствия теоретического закона распределения и экспериментально полученного 0,5.

На рисунке 3.3 изображен график экспоненты и гистограмма. Видно, что они почти согласуются. Это же подтверждается расчетом ч2.



Рисунок 3.3 - Сравнение теоретического закона распределения и экспериментально полученного

4 ВЫПОЛНЕНИЕ МОДЕЛИРОВАНИЯ РАБОТЫ СИСТЕМЫ ПРИ ВХОДНОМ СИГНАЛЕ N(t)

Для того чтобы получить выходной спектр нужно передаточную функцию системы Wобщ(p) (найденную в пункте 2) умножить на спектр входного сигнала (шума) N(p). Спектр сигнала n(t), находится с помощью прямого преобразования Фурье.

,(4.1)

Далее находятся спектры: амплитудный, фазовый, мнимой и реальной составляющих (см. рис. 4.1 - 4.4).

Рисунок 4.1 - Амплитудный спектр выходного шума

Рисунок 4.2 - Фазовый спектр выходного шума



Рисунок 4.3 - Реальная составляющая выходного шума

Рисунок 4.4 - Мнимая составляющая выходного шума

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



Рисунок 4.5 - Шум на выходе системы

Для того, чтобы получить автокорреляционную функцию системы (АКФ), используется теорема Винера-Хинчина. Эта теорема связывает между собой квадрат спектральной плотности сигнала и его корреляционную функцию через преобразование Фурье. По спектру выходного шума вычисляется энергетический спектр (формально это умножение на комплексно сопряженный). А после применения обратного преобразования Фурье получаем АКФ (см. рис. 4.6). Для построения ее графика было взято N/4 отсчетов.

Рисунок 4.6 - Автокорреляционная функция

Построение гистограммы выходного шума (см. рис. 4.7) проводится с помощью подпрограммы Rhist. Взято, как и в предыдущем задании 15 интервалов, для которых рассчитывается экспериментальное попадание в интервал с использованием подпрограммы для подсчета их площади по формуле Симпсона.

Рисунок 4.7 - Гистограмма выходного шума

Видно (из рис. 4.7), что закон распределения выходного шума уже не экспоненциальный. Для того, чтобы определить, к какому закону распределения ближе полученная гистограмма, нужно определить эксцесс и асимметрию. Асимметрия показывает насколько симметричен закон распределения

,(4.2)

где ni - отчеты выходного сигнала;

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

у - среднеквадратическое отклонение

,(4.3)

где DX - дисперсия.

Эксцесс показывает, на сколько островершинный полученный закон распределения и рассчитывается как

.(4.4)

Результаты расчета: mn=-2,52·10-5, у = 3,64 · 10-2, Ь = 0,32, е = 2,32. После нанесения этих точек на плоскость моментов (см. приложение А), видно, что ближайшими к полученному закону распределения являются распределения Релея и нормальный закон.

Плотность распределения вероятностей нормального закона имеет следующий вид

,(4.5)

Закон распределения Релея имеет следующий вид

,(4.6)

где х > 0.

Расчет теоретического количества попаданий отчетов случайного процесса в интервал гистограммы проводилось с помощью подпрограммы подсчета площади по формуле Симпсона от плотности распределения нормального закона распределения и распределения Рэлея.

Однако после проверки этих законов на критерий Пирсона, получаем:

а) для распределения Рэлея ч2 = 3601;

б) для нормального распределения ч2 = 1385.

Эти значения дают вероятность соответствия меньше 0,05.

На рис. 4.8 изображено совмещение теоретического количества попаданий в интервал гистограммы для распределения Рэлея и гистограммы выходного шума, а на рис. 4.9 совмещено теоретического количества попаданий в интервал гистограммы для нормального закона распределения и гистограммы выходного шума

Рисунок 4.8 - Совмещение теоретического количества попаданий в интервал гистограммы для распределения Рэлея и гистограммы выходного шума

Рисунок 4.9 - Совмещение теоретического количества попаданий в интервал гистограммы для нормального распределения и гистограммы выходного шума

ВЫВОДЫ

Заданный входной сигнал состоит из шести синусоид. В данном сигнале имеются по шесть мнимых и реальных отрицательных составляющих на частотах 75 Гц, 150 Гц, 55 Гц, 95 Гц, 130 Гц, 170 Гц. Однако реальных составляющих быть не должно. Их появление объясняется погрешностью вычисления. К тому же их значения достаточно малы - наибольшее из них -1,963 •10-3. В амплитудном спектре шесть составляющих, как это и должно быть, а их амплитуды соответствуют коэффициентам при синусоидах в сигнале. В фазовом спектре все составляющие компоненты с одинаковой фазой р/2.

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

Заданный шум распределен по экспоненциальному закону. При построении совместного графика теоретического распределения и гистограммы видно, что они хорошо согласуются. Это же подтверждается расчетами: ч2 = 11,68 и вероятность соответствия теоретического закона распределения и экспериментально полученного 0,5.

При прохождении шума через систему его форма искажается: форма его распределения уже не похожа на экспоненциальную. Для того, чтобы определить на какой закон распределения похож полученный на выходе шум, рассчитываются коэффициенты эксцесса и асимметрии. Далее эти точки наносятся на плоскость моментов, и определяется, какой из законов распределения ближайший. Наиболее подходящим оказался закон распределения Релея, но при проверке данного распределения по критерию Пирсона оказалось, что вероятность соответствия между экспериментально полученным законом распределения и распределением Релея меньше 0,05. это связано с тем, что Распределение Релея начинается с нуля, в отличии от экспериментально полученного. Поэтому также выполнилась оценка по критерию Пирсона на соответствие с нормальным законом распределения. В этом случае ч2 = 1385, что тоже дает вероятность соответствия меньше 0,05.

ПРИЛОЖЕНИЕ А

ПРИЛОЖЕНИЕ Б

!ПЕРВОЕ ЗАДАНИЕ - МОДЕЛИРОВАНИЕ СИГНАЛА S(t)=(1+а*соs2*pi*w*t)*(sinpi*f0*t+sin2*pi*f0*t),

!ПОСТРОЕНИЕ ГРАФИКА S(t) И СПЕКТРОВ СИГНАЛОВ.

INTEGER,PARAMETER:: N=1024, M=512, k=11

REAL dt/1.9531e-4/,df/5.0/

REAL S(N), t(N)

REAL SPECTR_A(M), SPECTR_Ph(M), SPECTR_Re(M), SPECTR_Im(M),f(M)

COMPLEX cS(N)

REAL a/0.75/, f0/150.0/, w/20.0/

REAL pi/3.1415/

REAL R1/5.1/, R2/0.1/, C1/68.0e-6/, L1/1.0e-2/

REAL R_1/10.0/, R_2/2.0/, L_1/5.0e-2/

REAL R__1/1000.0/, R__2/1.0/, C__1/47.0e-6/,L__1/1.0e-2/

REAL W9(M),W11(M),W18(M)

REAL FT(K)/0.0,260.0,520.0,780.0, 1040.0,1300.0,1560.0,1820.0,2080.0,2340.0,2600.0/

REAL AChH(K)/0.0, -6.54, -11.5, -15.8, -18.3, -19.9, -21.7, -24.0, -25.2, -28.2, -32.1/

REAL FChH(K)/0.0, -58.0, -72.7, -78.3, -81.2, -83.0, -84.2, -85.0, -85.7, -86.2, -86.6/

REAL AChHT(M), FChHT(M), AChHT_EDENICY(M), FChHT_RADIANY(M), WTT(M), W_USTROJSTVA(M)

REAL SpA2(M),SpPh2(M),SpRe2(M),SpIm2(M) ,SY(N)

INTEGER, PARAMETER::Nhist=15,NPOINT=50*NHIST+1,Ng=2

REAL*4 R, XSTART/0.37290981/,XKOEF/37.0/, E/2.7/,X(N),H(NHIST)

REAL MX/0.3/, DX, DX0, SIGMAX, XMAX, XMIN, TT(NHIST),FX(NPOINT),LYAMDA/3.3333/,XX(NPOINT)

REAL PLOSHA(NHIST),XF(120),XXX(NHIST),PLOSHA1(NHIST),XXXX(120),INTERVAL,DVA(NHIST,NG),TT1(NHIST,NG),XIXI

REAL XMAX2, XMIN2,TT2(NHIST),SxvyhR(N),RAKF(N),C,RAKF1(N),MXSh,DXSh, ALFAX, EX,SIGMAXSh,MX0Sh,DX0Sh,ALFAX0,EX0,H1(NHIST)

REAL XNORM(NPOINT), BNORM,ANORM,PLOTN(NPOINT),INTERVAL2,POPADANIE1(NHIST),DVA2(NHIST,NG),POPADANIE(NHIST),XNORMAL(NHIST,NG)

REAL XREL(NPOINT), BREL,AREL,PLOTN2(NPOINT),POPADANIE2(NHIST),DVA3(NHIST,NG),POPADANIE3(NHIST),XRELEJ(NHIST,NG),INTERVAL3

COMPLEX CP(N), CW9(N),CW11(N),CW18(N), CWT(M)

COMPLEX CW1(M), CW2(M), CW3(N), CW(N), CSY(N),CX(N),Sxvyh(N),SPECTRshuma_Re(M),SPECTRshuma_Im(M),SPECTRshuma_A(M),SPECTRshuma_Ph(M)

COMPLEX CJ/(0.0,1.0)/,AKF(N)

DO I=1,N

t(I)=REAL(I-1)*dt

S(I)=(1+a*COS(2*pi*w*t(I)))*(SIN(pi*f0*t(I))+SIN(2*pi*f0*t(I)))

END DO

!CALL DRAW(N,1,t,S,0,'LINE','МОДЕЛИРУЕМЫЙ ВХОДНОЙ СИГНАЛ')

!ACCEPT*

DO I=1,N

cS(I)=CMPLX(S(I),0.0)

END DO

CALL FFT (cS,N,0)

DO I=1,M

SPECTR_Re(I)=REAL(cS(I))

SPECTR_Im(I)=AIMAG(cS(I))

f(I)=REAL(I-1)*df

END DO

!CALL DRAW (M,1,f,SPECTR_Re,0,'STICK','РЕАЛЬНАЯ СОСТАВЛЯЮЩАЯ ВХОДНОГО СИГНАЛА')

!ACCEPT*

!CALL DRAW (M,1,f,SPECTR_Im,0,'STICK','МНИМАЯ СОСТАВЛЯЮЩАЯ ВХОДНОГО СИГНАЛА')

!ACCEPT*

DO I=1,M

SPECTR_A(I)=CABS(cS(I))

IF(ABS(REAL(cS(I)))<1.0e-3.AND.ABS(AIMAG(cS(I)))<1.0e-3) THEN

SPECTR_Ph(I)=0.0

ELSE

SPECTR_Ph(I)=ATAN2(AIMAG(cS(I)),REAL(CS(I)))

ENDIF

END DO

!CALL DRAW (M,1,f,SPECTR_A,0,'STICK','СПЕКТР ВХОДНОГО СИГНАЛА')

!ACCEPT*

!CALL DRAW (M,1,f,SPECTR_Ph,0,'STICK','ФАЗА ВХОДНОГО СИГНАЛА')

!ACCEPT*

! ВТОРОЕ ЗАДАНИЕ - ПОСТРОЕНИЕ МАТ. МОДЕЛИ СИСТЕМЫ.

DO I=1,M

CP(I)=2.0*pi*f(I)*CJ

ENDDO

!СОЗДАНИЕ МАССИВА ПЕРЕДАТОЧНЫХ ХАРАКТЕРИСТИК БЛОКОВ 9, 11, 18

DO I=1,M

CW9(I)=(CP(I)*L1+R2)/(CP(I)**2*R1*C1*L1+CP(I)*(R1*R2*C1+L1)+R1+R2)

CW11(I)=(CP(I)*L_1+R_2)/(CP(I)*L_1+R_1+R_2)

CW18(I)=(CP(I)**2*R__1*C__1*L__1+CP(I)*(R__1*R__2*C__1+L__1)+R__2)/(CP(I)**2*R__1*C__1*L__1+CP(I)*(R__1*R__2*C__1+L__1)+R__1+R__2)

END DO

!ПОЛУЧЕНИЕ АМПЛИТУДНых ХАРАКТЕРИСТИК БЛОКОВ

DO I=1,M

W9(I)=CABS(CW9(I))

W11(I)= CABS(CW11(I))

W18(I)= CABS(CW18(I))

END DO

! ПОСТРОЕНИЕ ПЕРЕДАТОЧНЫХ ХАРАКТЕРИСТИК БЛОКОВ И ЧХ Т-БЛОКА

!CALL DRAW (M,1,f,W9,0,'LINE','ПЕРЕДАТОЧНАЯ ХАРАКТЕРИСТИКА БЛОКА 9')

!ACCEPT*

!CALL DRAW (M,1,f,W11,0,'LINE','ПЕРЕДАТОЧНАЯ ХАРАКТЕРИСТИКА БЛОКА 11')

!ACCEPT*

!CALL DRAW (M,1,f,W18,0,'LINE','ПЕРЕДАТОЧНАЯ ХАРАКТЕРИСТИКА БЛОКА 18')

!ACCEPT*

!CALL DRAW (K,1,FT,AChH,0,'LINE','AЧХ Т БЛОКА')

!ACCEPT*

!CALL DRAW (K,1,FT,FChH,0,'LINE','ФЧХ Т БЛОКА')

!ACCEPT*

!!!!АПРОКСИМАЦИЯ!!!!

DO I=1,M

AChHT(I)=RINTERPOL(K,FT,AChH,F(I))

FChHT(I)=RINTERPOL(K,FT,FChH,F(I))

ENDDO

DO I=1,M

AChHT(I)=(AChHT(I))/20.0

AChHT_EDENICY(I)=10.0**AChHT(I)

FChHT_RADIANY(I)=(FChHT(I))*pi/180.0

ENDDO

! ПОСТРОЕНИЕ АПРОКСИМИРОВАНЫХ ЧХ Т-БЛОКА

!CALL DRAW (M,1,F,AChHT_EDENICY,0,'LINE','АПРОКСИМИРОВАНАЯ AЧХ Т-БЛОКА')

!ACCEPT*

!CALL DRAW (M,1,F,FChHT_RADIANY,0,'LINE','АПРОКСИМИРОВАНАЯ ФЧХ Т-БЛОКА')

!ACCEPT*

! СОЗДАНИЕ КОМПЛЕКСНОГО КОЕФИЦИЕНТА ПЕРЕДАЧИ Т-БЛОКА

DO I=1,M

CWT(I)=AChHT_EDENICY(I)*EXP(CJ*FChHT_RADIANY(I))

ENDDO

DO I=1,M

WTT(I)=CABS(CWT(I))

ENDDO

!CОЗДАНИЕ КОЕФИЦИЕНТА ПЕРЕДАЧИ ВСЕЙ СИСТЕМЫ

! БЛОК 11 И Т-БЛОК ПОСЛЕДОВАТЕЛЬНО

DO I=1,M

CW1(I)=CW11(I)*CWT(I)

ENDDO

! ПОЛУЧЕНЫЙ БЛОК 1 И БЛОК 18 С ОТРИЦАТЕЛЬНОЙ ОБР. СВ.

DO I=1,M

CW2(I)=CW1(I)/(1+CW1(I)*CW18(I))

ENDDO

! ПОЛУЧЕНЫЙ БЛОК 2 И БЛОК 9 ПОСЛЕДОВАТЕЛЬНО

DO I=1,M

CW3(I)=CW2(I)*CW9(I)

ENDDO

!МОДЕЛИРОВАНИЕ АЧХ ВСЕГО БЛОКА

DO I=1,M

W_USTROJSTVA(I)=CABS(CW3(I))

ENDDO

!CALL DRAW(M,1,F,W_USTROJSTVA,0,'LINES','ОБЩАЯ АМПЛИТУДО-ЧАСТОТНАЯ ХАРАКТЕРИСТИКА ВСЕГО БЛОКА')

!ACCEPT*

! СОЗДАНИЕ СИГНАЛА НА ВЫХОДЕ СИСТЕМЫ

DO I=1,M-1

CW3(M+I+1)=CONJG(CW3(M-I+1))

ENDDO

DO I=1,N

CSY(I)=CS(I)*CW3(I)

ENDDO

CALL FFT(CSY,N,1)

DO I=1,N

SY(I)=REAL(CSY(I))

ENDDO

CALL FFT(CSY,N,0)

DO I=1,M

SpRe2(I)=REAL(CSY(I))

SpIm2(I)=AIMAG(CSY(I))

ENDDO

!CALL DRAW(N,1,t,SY,0,'LINES','СИГНАЛ НА ВЫХОДЕ')

!ACCEPT*

!CALL DRAW(M,1,F,SpRe2,0,'STIKC','РЕАЛЬНАЯ СОСТАВЛЯЮЩАЯ СИГНАЛА НА ВЫХОДЕ')

!ACCEPT*

!CALL DRAW(M,1,F,SpIm2,0,'STIKC','МНИМАЯ СОСТАВЛЯЮЩАЯ СИГНАЛА НА ВЫХОДЕ')

!ACCEPT*

!СОЗДАНИЕ МАССИВОВ АМПЛИТУДНОЙ И ФАЗОВОЙ ХАРАКТЕРИСТИК

DO I=1,M

SpA2(I)=CABS(CSY(I))

IF((ABS(REAL(CSY(I)))<1E-3).and.ABS(IMAG(CSY(I)))<1E-3) THEN

SpPh2(I)=0.0

ELSE

SpPh2(I)=ATAN2(AIMAG(CSY(I)),REAL(CSY(I)))

ENDIF

ENDDO

!ПОСТРОЕНИЕ ЧХ ВЫХОДНОГО СИГНАЛА

!CALL DRAW(M,1,F,SpA2,0,'STIKC','АМПЛИТУДНО-ЧАСТОТНАЯ ХАРАКТЕРИСТИКА ВЫХОДНОГО СИГНАЛА')

!ACCEPT*

!CALL DRAW(M,1,F,SpPh2,0,'STIKC','ФАЗО-ЧАСТОТНАЯ ХАРАКТЕРИСТИКА ВЫХОДНОГО СИГНАЛА')

!ACCEPT*

! ТРЕТЬЕ ЗАДАНИЕ

!применение мультипликативного метода

DO I=1,N

XSTART=XSTART*XKOEF

XSTART=XSTART-AINT(XSTART)

R=XSTART

X(I)=-(LOG(1-R))/LYAMDA

! X(I) ПОЛУЧЕНО МЕТОДОМ ОБРАТНОЙ ФУНКЦИИ

ENDDO

!CALL DRAW(N,1,T,X,0,'LINE','ШУМ')

!ACCEPT*

! НАХОЖДЕНИЕ ДИСПЕРСИИ И СР. КВ. ОТКЛОНЕНИЯ

DX0=0.0

DO I=1,N

DX0=DX0+(X(I)-MX)**2

ENDDO

DX=DX0/REAL(N-1)

SIGMAX=SQRT(DX)

!PRINT *,'DX=',DX,'SKO=',SIGMAX

ACCEPT*

! НАХОЖДЕНИЕ МИН. И МАКС. ЭЛЕМЕНТОВ МАССИВА

XMIN=0.0

XMAX=X(1)

DO I=1,N

IF(X(I)>XMAX) XMAX=X(I)

ENDDO

!PRINT *,'XMIN=',XMIN, 'XMAX=',XMAX

! ПОСТРОЕНИЕ ГИСТОГРАМЫ

CALL RHIST(H,NHIST,X,N,XMIN,XMAX)

DO I=1,NHIST

TT(I)=XMIN+REAL(I-1)*(XMAX-XMIN)/(NHIST-1)

ENDDO

!ACCEPT*

!PRINT *,'TT=',TT

!ACCEPT*

!CALL DRAW(NHIST,1,TT,H,0,'BOX','ГИСТОГРАММА')

!ACCEPT*

!ПОСТРОЕНИЕ ФУНКЦИИ ЭКСПОНЕНТЫ

DO I=1,NPOINT

XX(I)=XMIN+((XMAX-XMIN)/(NPOINT-1))*REAL(I-1)

ENDDO

DO I=1,NPOINT

FX(I)=LYAMDA*EXP(-LYAMDA*XX(I))

ENDDO

!CALL DRAW(NPOINT,1,XX,FX,0,'LINES','ЭКСПОНЕНТА')

!ACCEPT*

INTERVAL=(XMAX-XMIN)/REAL(NPOINT)

CALL PL(PLOSHA,NHIST,INTERVAL,FX,NPOINT)

DO I=1,NHIST

PLOSHA1(I)=PLOSHA(I)*REAL(N)

ENDDO

SP=0.0

DO I=1,NHIST

SP=SP+PLOSHA(I)

ENDDO

PRINT *,PLOSHA,'SP',SP

ACCEPT*

!CALL DRAW (NHIST,1,TT,PLOSHA1,0,'LINES','ПЛОЩАДЬ')

!ACCEPT*

DO I=1,NHIST

DVA(I,2)=H(I)

DVA(I,1)=PLOSHA1(I)

TT1(I,1)=TT(I)

TT1(I,2)=TT(I)

ENDDO

!CALL DRAW(NHIST, NG, TT1, DVA, (/0,0/),(/'LINES','BOX'/),'СРАВНЕНИЕ')

!ACCEPT*

XIXI=0.0

DO I=1,NHIST

XIXI=XIXI+(PLOSHA1(I)-H(I))**2/PLOSHA1(I)

ENDDO

!PRINT *,'XIXI=', XIXI

!! 4 ЗАДАНИЕ

DO I=1,M

CX(I)=CMPLX(X(I)-MX,0.0)

ENDDO

CALL FFT (CX,N,0)

DO I=1,N

Sxvyh(I)=CX(I)*CW3(I)

ENDDO

!CALL DRAW (M,1,F,Sxvyh,0,'LINES','CПЕКТР ШУМА НА ВЫХОДЕ')

!ACCEPT*

DO I=1,M

SPECTRshuma_Re(I)=REAL(Sxvyh(I))

SPECTRshuma_Im(I)=AIMAG(Sxvyh(I))

END DO

!CALL DRAW (M,1,f,SPECTRshuma_Re,0,'STICK','РЕАЛЬНАЯ СОСТАВЛЯЮЩАЯ ВЫХОДНОГО ШУМА')

!ACCEPT*

!CALL DRAW (M,1,f,SPECTRshuma_Im,0,'STICK','МНИМАЯ СОСТАВЛЯЮЩАЯ ВЫХОДНОГО ШУМА')

!ACCEPT*

DO I=1,M

SPECTRshuma_A(I)=CABS(Sxvyh(I))

IF(ABS(REAL(Sxvyh(I)))<1.0e-3.AND.ABS(AIMAG(Sxvyh(I)))<1.0e-3) THEN

SPECTRshuma_Ph(I)=0.0

ELSE

SPECTRshuma_Ph(I)=ATAN2(AIMAG(Sxvyh(I)),REAL(Sxvyh(I)))

ENDIF

END DO

!CALL DRAW (M,1,f,SPECTRshuma_A,0,'STICK','АМПЛИТУДНЫЙ СПЕКТР ВЫХОДНОГО ШУМА')

!ACCEPT*

!CALL DRAW (M,1,f,SPECTRshuma_Ph,0,'STICK','ФАЗА ВЫХОДНОГО ШУМА')

!ACCEPT*

DO I=1,N

AKF(I)=Sxvyh(I)*CONJG(Sxvyh(I))

ENDDO

CALL FFT(AKF,N,1)

DO I=1,N

RAKF(I)=REAL(AKF(I))

ENDDO

C=AKF(1)

DO I=1,N

RAKF1(I)=RAKF(I)/C

ENDDO

!CALL DRAW (N/4,1,T,RAKF1,0,'LINE','AKF')

!ACCEPT*

CALL FFT(Sxvyh,N,1)

DO I=1,N

SxvyhR(I)=REAL(Sxvyh(I))

ENDDO

XMIN2=SxvyhR(1)

XMAX2=SxvyhR(1)

DO I=1,N

IF(SxvyhR(I)>XMAX2) XMAX2=SxvyhR(I)

IF(SxvyhR(I)<XMIN2) XMIN2=SxvyhR(I)

ENDDO

PRINT *,'XMIN2=',XMIN2, 'XMAX2=',XMAX2

CALL RHIST(H1,NHIST,SxvyhR,N,XMIN2,XMAX2)

DO I=1,NHIST

TT2(I)=XMIN2+REAL(I-1)*(XMAX2-XMIN2)/(NHIST-1)

ENDDO

!CALL DRAW (M,1,T,SxvyhR,0,'LINE','ШУМ НА ВЫХОДЕ СИСТЕМЫ')

!ACCEPT*

!CALL DRAW(NHIST,1,TT2,H1,0,'BOX','ГИСТОГРАММА ВЫХОДНОГО ШУМА')

!ACCEPT*

!MX0Sh=0.0

DO I=1,N

MX0Sh=MX0Sh+SxvyhR(I)

ENDDO

MXSh=MX0Sh/REAL(N)

DX0Sh=0.0

DO I=1,N

DX0Sh=DX0Sh+(SxvyhR(I)-MXSh)**2

ENDDO

DXSh=DX0Sh/REAL(N-1)

SIGMAXSh=SQRT(DXSh)

ALFAX0=0.0

DO I=1,N

ALFAX0=ALFAX0+((SxvyhR(I)-MXSh)**3)

ENDDO

ALFAX=ALFAX0/REAL(N*SIGMAXSh**3)

EX0=0.0

DO I=1,N

EX0=EX0+((SxvyhR(I)-MXSh)**4)

ENDDO

EX=EX0/(REAL(N)*(SIGMAXSh**4))-3

PRINT*, 'MXSh=',MXSh, ' DXSh=',DXSh, 'ALFAX=', ALFAX, 'EX=', EX, 'SIGMAXSh=', SIGMAXSh

DO I=1,NPOINT

! ПОСТРОЕНИЕ ПЛОТНОСТИ НОРМАЛЬНОГО РАССПРЕДЕЛЕНИЯ

XNORM(I)=XMIN2+REAL(I-1)*(XMAX2-XMIN2)/REAL(NPOINT-1)

ENDDO

BNORM=MXSh+3.0*SIGMAXSh

ANORM=MXsH-3.0*SIGMAXSh

DO I=1,NPOINT

PLOTN(I)=(1.0/(SQRT(2.0*PI)*SIGMAXSh))*EXP(-(1.0/(2.0*SIGMAXSh**2))*(XNORM(I)-MXSh)**2)

ENDDO

!CALL DRAW(NPOINT,1,XNORM,PLOTN,0,'LINE','НОРМАЛЬНОЕ РАСПРЕДЕЛЕНИЕ')

!ACCEPT*

! ФОРМИРОВАНИЕ МАССИВА КОЛИЧЕСТВА ПОПАДАНИЙ В ИНТЕРВАЛ ГИСТОГРАММЫ

INTERVAL2=(XMAX2-XMIN2)/REAL(NPOINT)

CALL PL(POPADANIE,NHIST,INTERVAL2,PLOTN,NPOINT)

DO I=1,NHIST

POPADANIE1(I)=POPADANIE(I)*REAL(N)

ENDDO

DO I=1,NHIST

DVA2(I,2)=H1(I)

DVA2(I,1)=POPADANIE1(I)

XNORMAL(I,1)=TT2(I)

XNORMAL(I,2)=TT2(I)

ENDDO

!CALL DRAW(NHIST, NG, XNORMAL, DVA2, (/0,0/),(/'LINES','BOX'/),'СРАВНЕНИЕ')

!ACCEPT*

XIXINORM=0.0

DO I=1,NHIST

XIXINORM=XIXINORM+(POPADANIE1(I)-H1(I))**2/POPADANIE1(I)

ENDDO

PRINT *,'XIXINORM=', XIXINORM

! ПОСТРОЕНИЕ ПЛОТНОСТИ РАССПРЕДЕЛЕНИЯ РЕЛЛЕЯ

BREL=XMAX2

AREL=0.0

DO I=1,NPOINT

XREL(I)=AREL+REAL(I-1)*(BREL-AREL)/REAL(NPOINT-1)

ENDDO

DO I=1,NPOINT

PLOTN2(I)=(XREL(I)/SIGMAXSh**2)*EXP(-XREL(I)**2/(2.0*SIGMAXSh**2))

ENDDO

!CALL DRAW(NPOINT,1,XREL,PLOTN2,0,'LINE','РАСПРЕДЕЛЕНИЕ РЕЛЕЯ')

!ACCEPT*

!ФОРМИРОВАНИЕ МАССИВА КОЛИЧЕСТВА ПОПАДАНИЙ В ИНТЕРВАЛ ГИСТОГРАММЫ

INTERVAL3=(BREL-AREL)/REAL(NPOINT)

CALL PL(POPADANIE2,NHIST,INTERVAL3,PLOTN2,NPOINT)

DO I=1,NHIST

POPADANIE3(I)=POPADANIE2(I)*REAL(N)

ENDDO

!CALL DRAW(NHIST,1,TT2,POPADANIE2,0,'LINE','ПЛОЩАДЬ')

!ACCEPT*

do I=1,NHIST

DVA3(I,2)=H1(I)

DVA3(I,1)=POPADANIE3(I)

XRELEJ(I,1)=TT2(I)

XRELEJ(I,2)=TT2(I)

ENDDO

!CALL DRAW(NHIST, NG, XREL, DVA3, (/0,0/),(/'LINES','BOX'/),'СРАВНЕНИЕ')

!ACCEPT*

XIXIREL=0.0

DO I=1,NHIST

XIXIREL=XIXIREL+(POPADANIE3(I)-H1(I))**2/POPADANIE3(I)

ENDDO

PRINT *,'XIXIREL=', XIXIREL

END

ПОДПРОГРАММА ИНТЕРПОЛЯЦИИ

REAL FUNCTION RINTERPOL(L, X_ZAD, f_ZAD, X)

INTEGER L

REAL X_ZAD(L), f_ZAD(L), X

REAL CH(L),ZN(L),WT

DO I=1,L

CH(I)=1.0

ZN(I)=1.0

DO J=1,L

IF(I/=J) CH(I)=CH(I)*(X-X_ZAD(J))

IF(I/=J) THEN

ZN(I)=ZN(I)*(X_ZAD(I)-X_ZAD(J))

ENDIF

ENDDO

ENDDO

WT=0.0

DO I=1,L

WT=WT+f_ZAD(I)*CH(I)/ZN(I)

ENDDO

RINTERPOL=WT

RETURN

END FUNCTION RINTERPOL

ПОДПРОГРАММА ПОДСЧЕТА ПЛОЩАДИ

SUBROUTINE PL(S,N1,H,Y,N)

INTEGER N,N1

REAL Y(N)

REAL SUMMA,H,S(N1)

DO J=2,N-1,2

Y(J)=4.0*Y(J)

ENDDO

DO I=1,N,2

Y(I)=2.0*Y(I)

ENDDO

DO J=1,N1

SUMMA=0.0

DO I=(J-1)*INT(N/N1)+1,J*INT(N/N1)+1

Y((J-1)*INT(N/N1)+1)=Y((J-1)*INT(N/N1)+1)/2.0

Y(J*INT(N/N1)+1)=Y(J*INT(N/N1)+1)/2.0

SUMMA=SUMMA+Y(I)

ENDDO

S(J)=SUMMA*H/3.0

ENDDO

RETURN

END

Размещено на Allbest.ru