SM10
Программа SM10 вычисляет сумму бесконечного ряда с помощью преобразования Эйлера, которое заключаетсмя в следующем. Если ряд в левой части равенства
∞∞
∑ ai = ∑ (∆ia0)/2i+1
i=0i=0
сходится, то это равенство справедливо и ряд в правой части равенства сходится быстрее. Здесь ∆ia0 обозначает разность вперед i-го порядка для a0
∆a0 = a1 - a0
∆1a0 = ∆a1 - ∆a0
∆2a0 = ∆1a1 - ∆1a0
. . . . . . . . . .
Программа особенно эффективна в случае медленно сходящихся, а также расходящихся знакопеременных рядов. Суммирование прекращается, как только Tim раз подряд абсолютные значения членов преобразованного ряда будут меньше, чем Eps.
Вызов программы
SM10(Fct, Eps, Tim)
Параметры программы
Fct, Eps, Tim - входные параметры;
Real SM10 - возвращаемое значение;
Real Fct(Integer i) - функция, вычисляющая значения членов ряда;
Real Eps - погрешность вычисления суммы;
Integer Tim - параметр повторения.
Пример
!Вычисление суммы знакопеременного ряда program TestSM10 use NML implicit none real, parameter:: pi=3.1415926, Eps=1.E-8 integer, parameter:: Tim=4 real:: S, E !begin S=SM10(Pi4, Eps, Tim) E=0.25*pi print 10, S, E 10 format(/' Summa =',E18.10/' pi/4 =',E18.10) !end contains real function Pi4(i) integer, intent(in):: i Pi4=1.0/float(2*i+1) if(mod(i,2) /= 0) Pi4=-Pi4 return end function Pi4 end program TestSM10 Summa = 0.7853981256E+00 pi/4 = 0.7853981256E+00
DFT
Программа DFT вычисляет дискретное (прямое и обратное) быстрое преобразование Фурье с прореживанием по времени. Алгоритм Кули-Тьюки.
Прямое и обратное дискретное преобразование Фурье задаётся формулами
где s(n) - исходный сигнал, S(k) - его Фурье-образ, j - мнимая единица. Предполагается, что сигнал s(n) - периодический, в этом случае Фурье-образ будет тоже периодическим с тем же периодом.
Основная идея быстрого преобразования Фурье (БПФ) состоит в разбиении отсчётов сигнала на две части и вычислении БПФ уже для каждой из частей меньшей размерности N/2 (прореживание по времени). Во втором варианте БПФ (прореживание по частоте) на две части разбивается спектр сигнала. Подробное описание обоих алгоритмов можно найти в многочисленной литературе.
[Б14]
Вызов программы
call DFT(N, S, Error, Dir)
Параметры программы
N, Dir - входные параметры;
S - входной и выходной параметр;
Error - выходной параметр;
Integer N - количество отсчётов сигнала или Фурье-образа, N должно быть степенью числа 2, большей нуля;
Complex S[0:N-1] - массив с отсчётами;
Integer Dir - направление преобразования, Dir>0 или отсутствует - прямое преобразование Фурье, Dir<0 - обратное преобразование Фурье;
Integer Error - индикатор ошибки, Error=0 в случае нормального завершения программы, Error>0 если размерность массива S не является степенью числа 2.
Пример
! Прямое и обратное преобразование Фурье треугольного сигнала
program DigitalFourierTransform
use NML
implicit none
! Variables
integer, parameter:: n=8, n1=n-1
complex S(0:n1)
real Arg(0:n1), X(0:n1)
real a, b, h, t
integer k, err
data a/0.0/, b/2.0/
! Body of DigitalFourierTransform
h=(b-a)/n;
do k=0, n1
Arg(k)=a+h*k
if(k<=n/2) X(k)=h*k
if(k>n/2) X(k)=X(k-1)-h
S(k)=cmplx(X(k), 0.0)
enddo
print *, 'Original signal'
print 10, (k, Arg(k), X(k), k=0, n1)
call DFT(n, S, err)
print *, 'Error=',err
print *, 'The direct transform signal'
print 11, (k, Arg(k), Real(S(k)), Aimag(S(k)), k=0, n1)
call DFT(n, S, err, -1)
print *, 'The inverse transform signal'
print 11, (k, Arg(k), Real(S(k)), Aimag(S(k)), k=0, n1)
10 format(I2,' Arg=',F4.2,F10.5)
11 format(I2,' Arg=',F4.2,' Re=',F10.6,' Im=',F10.6)
end program DigitalFourierTransform
Original signal
0 Arg=0.00 0.00000
1 Arg=0.25 0.25000
2 Arg=0.50 0.50000
3 Arg=0.75 0.75000
4 Arg=1.00 1.00000
5 Arg=1.25 0.75000
6 Arg=1.50 0.50000
7 Arg=1.75 0.25000
Error= 0
The direct transform signal
0 Arg=0.00 Re= 4.000000 Im= 0.000000
1 Arg=0.25 Re= -1.707107 Im= 0.000000
2 Arg=0.50 Re= 0.000000 Im= 0.000000
3 Arg=0.75 Re= -0.292893 Im= -0.000000
4 Arg=1.00 Re= 0.000000 Im= 0.000000
5 Arg=1.25 Re= -0.292893 Im= 0.000000
6 Arg=1.50 Re= 0.000000 Im= 0.000000
7 Arg=1.75 Re= -1.707107 Im= 0.000000
The inverse transform signal
0 Arg=0.00 Re= 0.000000 Im= 0.000000
1 Arg=0.25 Re= 0.250000 Im= -0.000000
2 Arg=0.50 Re= 0.500000 Im= 0.000000
3 Arg=0.75 Re= 0.750000 Im= -0.000000
4 Arg=1.00 Re= 1.000000 Im= 0.000000
5 Arg=1.25 Re= 0.750000 Im= 0.000000
6 Arg=1.50 Re= 0.500000 Im= 0.000000
7 Arg=1.75 Re= 0.250000 Im= 0.000000
CNM
Программа CNM вычисляет биномиальные коэффициенты Cnm = n!/(m!(n-m)!) по рекуррентной формуле Cni+1 = (n-i)/(i+1)Cni для i = 0, 1, ..., m-1, причем Cn0 = 1.
Вызов программы
CNM(n, m)
Параметры программы
Integer n, m - входные параметры;
Integer CNM - возвращаемое значение.