Alg i metodi vichisl / Teorija / lek_num_met
.pdf
Вычисления всех 
и ci требуют O(N) действий, три преобразования Фурье требуют O(Nlog(N)) действий, перемножение результатов преобразований Фурье требует O(N) действий, вычисление всех bi зная значения свертки требует O(N) действий. Итого для дискретного преобразования Фурье требуется O(Nlog(N)) действий для любого N.
Этот алгоритм быстрого преобразования Фурье может работать над кольцом только когда известны первообразные корни из единицы степеней 2N и 2k.
2.4 Принцип работы Быстрого преобразования Фурье
БПФ работает с разложением N-точечного временного промежутка сигнала на N временных сигнальных промежутков, каждый из которых состоял из одной точки. Второй шаг заключается в расчете N частотных спектров, соответствующего этим временным сигнальным промежуткам. После этого, на основе N спектров синтезируется единый частотный спектр.
Рисунок 10.1 Процедура разбиения сигналапоказывает пример разделения временной области с использованием БПФ. В этом примере 16-ти точечный сигнал разлагается на четыре отдельные компоненты. Первый этап заключается в разбиения 16-ти точечного сигнала на два сигнала по 8 точек. На втором этапе происходит разложение данных сигналов на четыре сигнала по 4 точки. Данная процедура продолжается до тех пор, пока не будет сформировано N сигналов состоящих из одной точки. Чересстрочное разложения используется при разделении выборки сигнала на две выборки. Сигнал разделяется как при четном, так и при нечетном количестве отсчетов в выборке. В системе проводится 2N этапов в разложении, т.е. на 16–ти точечный сигнал (24) предусматривает 4 этапа, 512 точки (27) требует 7 этапов, 4096 точки (212) предусматривает 12 этапов и т.д. Запомним, это значение - оно будет много раз встречаться в этой главе.
Рисунок 10.1 Процедура разбиения сигнала
Теперь, когда понятна структура разложения, ее можно значительно упростить. Разложения является не более чем перераспределение отсчетов сигнала. Это иллюстрирует Рисунок 10.2 Процедура перераспределения отсчетов сигнала Слева представлены номера отсчетов первоначального сигнала, вместе с их двоичными эквивалентами. Справа, измененные номера, а также их двоичные эквиваленты. Идея состоит в том, что двоичные числа являются инверсией друг друга. Например, отсчёт №3 (0011) можно получить из отсчёта №12 (1100). Отсчёт № 14 (1110) получается из отсчета № 7 (0111), и так далее. Разбиение временной области в БПФ обычно, проводится с использование алгоритма сортировки инвертированных бит. Это предполагает изменение N временных областей, с подсчетом в двоичном виде и переворачиванием бит слева на право.
Рисунок 10.2 Процедура перераспределения отсчетов сигнала !!!
Следующим шагом алгоритма БПФ является нахождение частоты спектров одноточечного сигналов, определенного во временной области. Это делается элементарно: спектр одноточечного сигнала равен самому себе. Это означает, что на этом шаге ничего не надо. Хотя не стоит забывать, что каждый из одноточечных сигналов теперь представляет собой спектр частот, а не сигнал во временной области
Рисунок 10.5 Схема "Бабочки" !!!
Рисунок 10.6 Блок-схема программы, реализующая БПФпоказывает структуру БПФ. Разбиение во временной области выполняется с помощью ранее рассмотренного алгоритма сортировки.
Синтез в частотной области требует трех циклов. Внешний цикл проводится Log2N раз. Средний цикл движется по каждой отдельной частоте спектра. Внутренний цикл использует бабочку для расчета точек в каждой частоты спектра (т.е. циклы через образцы внутри одного окна на рис.1). Накладных клетки на рис.6 определяют начало и окончания индексов в циклах, а также расчета необходимых синусоиды бабочек. Блоксхема программы, реализующая БПФ приведена на Рисунок 10.6 Блок-схема программы, реализующая БПФ
Рисунок 10.6 Блок-схема программы, реализующая БПФ !!!
2.5 Программа, реализующая БПФ
Реализация БПФ на языке программирования C
/*************************************************************************
Алгоритм проводит быстрое преобразование Фурье комплексной функции, заданной nn отсчетами на действительной оси.
В зависимости от переданных параметров, может выполняться как прямое, так и обратное преобразование.
Входные параметры:
nn- Число значений функции. Должно быть степенью двойки. Алгоритм не проверяет правильность переданного значения.
a - array [0 .. 2*nn-1] of Real
Значения функции. I-ому значению соответствуют элементы a[2*I] (вещественная часть) и a[2*I+1] (мнимая часть).
InverseFFT
-направление преобразования.
True, если обратное, False, если прямое.
Выходные параметры:
a - результат преобразования.
*************************************************************************/
void fastfouriertransform(ap::real_1d_array& a, int nn, bool inversefft)
{
int ii; int jj; int n; int mmax; int m; int j; int istep; int i; int isign;
double wtemp; double wr; double wpr; double wpi; double wi; double theta; double tempr; double tempi;
if( inversefft )
{
isign = -1;
}
else
{
isign = 1;
}
n = 2*nn; j = 1;
for(ii = 1; ii <= nn; ii++)
{
i = 2*ii-1; if( j>i )
{
tempr = a(j-1); tempi = a(j); a(j-1) = a(i-1); a(j) = a(i);
a(i-1) = tempr; a(i) = tempi;
}
m = n/2; while(m>=2&&j>m)
{
j = j-m; m = m/2;
}
j = j+m;
}
mmax = 2; while(n>mmax)
{
istep = 2*mmax;
theta = 2*ap::pi()/(isign*mmax); wpr = -2.0*ap::sqr(sin(0.5*theta)); wpi = sin(theta);
wr = 1.0; wi = 0.0;
for(ii = 1; ii <= mmax/2; ii++)
{
m = 2*ii-1;
for(jj = 0; jj <= (n-m)/istep; jj++)
{
i = m+jj*istep; j = i+mmax;
tempr = wr*a(j-1)-wi*a(j); tempi = wr*a(j)+wi*a(j-1); a(j-1) = a(i-1)-tempr; a(j) = a(i)-tempi;
a(i-1) = a(i-1)+tempr; a(i) = a(i)+tempi;
}
wtemp = wr;
wr = wr*wpr-wi*wpi+wr; wi = wi*wpr+wtemp*wpi+wi;
}
mmax = istep;
}
if( inversefft )
{
for(i = 1; i <= 2*nn; i++)
{
a(i-1) = a(i-1)/nn;
}
}
}
Реализация БПФ на языке программирования Pascal
procedure LinearFDFT (Re, Im: PReal; Count: Word; Direct: ShortInt); var
I: Word;
K: Real;
PIm, PRe: PReal;
procedure LFDFT (SrcRe, SrcIm: PReal; Cnt: Word); var
EvenRe, OddRe, PEvenRe, POddRe, PRe, PSrcRe: PReal; …; //то же – для мнимых
Factor: Real;
HalfCnt, I, Size: Word;
HEvenIm, HEvenRe, HOddIm, HOddRe: THandle; X, Y, WIm, WRe: Real;
begin
PSrcRe := SrcRe; …; if Cnt = 2 then
begin //тривиальная операция «бабочка»
Inc (PSrcRe); …;
X := SrcRe^; Y := PSrcRe^;
SrcRe^ := X + Y; PSrcRe^ := X – Y; X := SrcIm^; Y := PSrcIm^;
SrcIm^ := X + Y; PSrcIm^ := X – Y; end else
begin //переупорядочение и рекурсивный вызов для полпоследовательностей
Factor := K / Cnt; HalfCnt := Cnt div 2;
Size := HalfCnt * SizeOfReal + 1;
HEvenRe := GlobalAlloc (GMem_Fixed, Size); …; HOddRe := GlobalAlloc (GMem_Fixed, Size); …; EvenRe := GlobalLock (HEvenRe); …;
OddRe := GlobalLock (HOddRe); …; PEvenRe := EvenRe; …;
POddRe := OddRe; …;
//Разбивка на (не)чётные подпоследовательности for I := 0 to Cnt – 1 do
begin
if Odd (I) then begin
POddRe^ := PSrcRe^; …; Inc (POddRe); …;
end else begin
PEvenRe^ := PSrcRe^; …; Inc (PEvenRe); …;
end;
Inc (PSrcRe); …; end;
//Рекурсивная обработка (не)чётных подпоследовательностей
LFDFT (EvenRe, EvenIm, HalfCnt);
LFDFT (OddRe, OddIm, HalfCnt);
//Сборка обработанных подпоследовательностей
PSrcRe := SrcRe; …; PRe := SrcRe; …;
Inc (pRe, HalfCnt); …; PEvenRe := EvenRe; …; POddRe := OddRe; …;
for I := 0 to HalfCnt – 1 do begin
WRe := Cos (Factor * I); WIm := – Sin (Factor * I); PSrcRe^ := POddRe^ * WRe – POddIm^ * WIm;
PSrcIm^ := POddRe^ * WIm + POddIm^ * WRe; PRe^ := pSrcRe^; …;
PSrcRe^ := pEvenRe^ + PSrcRe^; …; PRe^ := PEvenRe^ – PRe^; …;
Inc (PSrcRe); …; Inc (PRe); …; Inc (PEvenRe); …; Inc (POddRe); …;
end;
GlobalUnlock (HEvenRe); …; GlobalUnlock (HOddRe); …; GlobalFree (HEvenRe); …; GlobalFree (HOddRe); …;
end;
end; begin
K := 2 * Pi * Direct;
LFDFT (Re, Im, Count); if Direct < 0 then
begin //Дополнительное деление на n для обратного преобразования
PRe := Re; …;
for I := 1 to Count do begin
PRe^ := PRe^ / Count; …; Inc (PRe); …;
end;
end;
end;
end.
Глава 11 АЛГОРИТМЫ ГЕНЕРАЦИИ СЛУЧАЙНЫХ ЧИСЕЛ С РАЗЛИЧНЫМИ
ЗАКОНАМИ РАСПЕРЕДЕЛЕНИЯ
§1. Генерация случайных чисел с нормальным законом распределения.
1.1 Алгоритмы моделирования случайных чисел
Дискретные случайные величины.
Слова "случайная величина" в обыденном смысле употребляют тогда, когда хотят подчеркнуть, что неизвестно, каким будет конкретное значение этой величины. Причем иногда за этими словами скрывается просто незнание, какова эта величина.
Итак, чтобы задать случайную величину, надо указать, какие значения она может принимать, и каковы вероятности этих значений.
Определение
Случайная величина Х называется дискретной, если она может принимать дискретное множество значений х1, х2, …, хn.
Формально случайная дискретная величина Х определяется таблицей
x |
x |
2 |
x |
|
(11.1) |
1 |
|
|
n |
||
X = p |
p |
2 |
p |
|
|
1 |
|
|
n |
|
|
где х1, х2, …, хn - возможные значения величины Х; |
|
|
|
|
|
р1, р2, …, рn - соответствующие вероятности. |
|
|
|
|
|
Точнее говоря, вероятность P{X=хi } того, что случайная величина Х примет значение хi, равна: |
|
||||
P{X = xi }= pi |
(11.2) |
|
|||
!!!!Таблица (1.1) называется распределением случайной дискретной величины. Числа х1, х2, …, хn могут быть вообще говоря, любыми. Однако вероятности р1, р2, …, рn должны удовлетворять двум условиям:
pi |
> 0 |
(11.3) |
и |
|
|
p1 + p2 |
+… + pn = 1 |
(11.4) |
Последнее условие означает, что Х обязана в каждом случае принять одно из значений х1, х2, …, хn.
Кроме распределения случайной величины, которая является исчерпывающей характеристикой, вводятся числовые характеристики, основными среди которых являются математическое ожидание и дисперсия. Получение случайных величин на ЭВМ.
Сама постановка вопроса "получение случайных чисел на ЭВМ" иногда вызывает недоумение: ведь все, что делает компьютер, должно быть заранее запрограммировано; откуда же может появиться случайность? Вопрос о том, можно ли с их помощью описать какое-либо явление природы, решается опытным путем. Такое описание всегда является приближенным. Более того, случайная величина, которая вполне удовлетворительно описывает какую-то физическую величину в одном классе явлений, может оказаться плохой характеристикой этой же величины при исследовании других явлений.
Обычно различают три способа получения случайных величин:
•из заранее составленных таблиц случайных чисел;
•физические генераторы случайных чисел;
•с помощью формул (генераторов или датчиков) псевдослучайных чисел.
Поскольку "качество" используемых в имитационном моделировании случайных чисел проверяется с помощью специальных тестов, можно не интересоваться тем, как эти числа получены: лишь бы они удовлетворяли принятой системе тестов.
Числа, получаемые по какой-либо формуле и имитирующие значения случайной величины X, называются псевдослучайными числами. Под словом "имитирующие" подразумевается, что эти числа удовлетворяют ряду тестов так, как если бы они были значениями этой случайной величины.
Основой или «сырьем» для моделирования случайных величин с заданным законом распределения являются так называемые базовые случайные числа.
Определение. Совокупность {Ri}, i=1,2, … независимых равномерно распределенных на отрезке [0,1] случайных величин называется последовательностью базовых случайных чисел.
Мы называем эти числа псевдослучайными потому, что фактически они остаются полностью детерминированными в том смысле, что если каждое обращение к соответствующей формуле (точнее, к алгоритму) начинается с одними и теми же исходными данными (константами и начальными значениями), то на выходе получаются одинаковые последовательности чисел R.
В настоящее время почти все стандартные библиотечные программы вычисления равномерных случайных
чисел основаны на конгруэнтных1 методах, разработанных Лемером. |
|
Основная формула мультипликативного конгруэнтного метода Лемера имеет вид: |
|
Ri +1 = aR i (mod m ) |
(11.5) |
где а и m - неотрицательные целые числа.
Согласно этому выражению, нужно взять случайное число Ri, умножить его на постоянный коэффициент а и взять модуль полученного числа по m (т.е. разделить на аRi и остаток считать как Ri+1). Поэтому для вычисления (или генерирования) последовательности Ri нам необходимы начальные значения R0, множитель а и модуль m. Выбираются а, R0 и m так, чтобы обеспечить максимальную длину (или, как говорят период) неповторяющейся последовательности Ri и минимальную корреляцию между генерируемыми числами. Базовые случайные числа позволяют генерировать новые случайные последовательности, подчиняющиеся любому закону распределения.
Существует два основных пути преобразования базовых случайных чисел {Ri}, в случайные числа {yi}, распределенные по заданному закону распределения.
Один из них, который называется методом инверсии, состоит в реализации определенных арифметических операций над базовым числом Ri, чтобы получить уi.
Второй метод основывается на моделировании условий соответствующей предельной теоремы теории вероятностей. Кроме указанных двух основных подходов можно также выделить эвристические способы генерирования случайных чисел.
Моделирование нормальной случайной величины на основе центральной предельной теоремы.
Нормальное (или гауссово) распределение (рисунок 11.1)—это, несомненно, один из наиболее важных и часто используемых в имитационном моделировании видов непрерывных распределений.
Рисунок 11.1 Нормальная плотность распределения !!!
Плотность вероятности нормальной плотности записывается так:
p(x) = |
|
1 |
e− |
(x−Mi )2 |
|
σ |
2δ 2 |
, |
|||
|
2π |
|
|
|
|
где Мх – математическое ожидание; σ - среднеквадратическое отклонение.
Функция распределения нормальной случайной величины равна:
|
|
l |
x |
(x−M x )2 |
|
P(x) = |
|
∫e− |
|
dx |
|
|
2δ 2 |
||||
|
σ |
2π −∞ |
|
|
|
Поэтому алгоритмы моделирования нормальных случайных чисел базируются на предельных теоремах теории вероятностей.
Центральная предельная теорема говорит о том, что сумма n одинаково распределенных независимых случайных величин x со средним Мх и дисперсией Dx стремится к нормально распределенной величине с параметрами nМх и nDx при бесконечном увеличении n.
Следствием теоремы является, в частности, и то, что для получения нормальной выборки, можно воспользоваться базовыми случайными числами R.
Идея алгоритма состоит в следующем. Определим новую случайную величину s в виде суммы базовых чисел
Ri, ( i=1,2,3, …,n):
