Известно:
- Код: Выделить всё
' Построение взаимного амплитудного спектра, функций разности фаз и
' взаимной когерентности двух случайных процессов по оценкам их
' автоспектров и взаимных спектров.
' Замечания по использованию:
' 1. с целью экономии числа используемых массивов результаты можно
' получать на месте исходной информации, а именно, допустимы
' совпадения параметров: DPH = SXY2, AV = SXY1, COH = SY.
Public Sub RSF1R(ByRef SX() As Double, _
ByRef SY() As Double, _
ByRef SXY1() As Double, _
ByRef SXY2() As Double, _
ByVal NC As Long, _
ByRef DPH() As Double, _
ByRef AV() As Double, _
ByRef COH() As Double)
' Входные параметры:
' SX - одномерный массив длины NC, содержащий заданные значения
' автоспектра процесса X (тип: вещественный);
' SY - одномерный массив длины NC, содержащий заданные значения
' автоспектра процесса Y (тип: вещественный);
' SXY1 - одномерный массив длины NC, содержащий заданные значения
' вещественной части взаимного спектра процессов X, Y
' (тип: вещественный);
' SXY2 - одномерный массив длины NC, содержащий заданные значения
' мнимой части взаимного спектра процессов X, Y (тип:
' вещественный);
' NC - число заданных значений автоспектров и взаимного спектра
' процессов X, Y (тип: целый);
' Выходные параметры:
' DPH - одномерный массив длины NC, содержащий вычисленные значения
' функции разности фаз процессов X, Y (тип: вещественный);
' AV - одномерный массив длины NC, содержащий вычисленные значения
' взаимного амплитудного спектра процессов X, Y (тип:
' вещественный);
' COH - одномерный массив длины NC, содержащий вычисленные значения
' функции взаимной когерентности процессов X, Y (тип:
' вещественный).
Dim I As Long
Dim Q As Double, R As Double, T As Double
For I = 1 To NC
Q = SXY1(I)
R = SXY2(I)
DPH(I) = Atn(R / Q)
T = Abs(Q): If Abs(R) > T Then T = Abs(R) ' T=AMAX1(ABS(Q),ABS(R))
Q = Q / T
R = R / T
Q = T * Sqr(Q * Q + R * R)
AV(I) = Q
COH(I) = Q * Q / (SX(I) * SY(I))
Next I
End Sub
' Построение автопериодограмм и взаимных периодограмм двух вещественных
' случайных процессов по одному отрезку их реализаций.
' Замечания по использованию:
' 1. После выхода из подпрограммы RSS1R значения исходных случайных процессов
' в массивах X1, Y1 не сохраняются.
' 2. B частном случае построения автопериодограммы одного случайного процесса
' X в массив Y1 следует заслать нули, при этом допустимы совпадения
' параметров: PX = PY = PXY1 = PXY2.
' 3. B частном случае построения только автопериодограмм двух случайных
' процессов X, Y без вычисления их взаимной периодограммы допустимы
' совпадения параметров: PXY1 = PXY2.
Public Sub RSS1R(ByRef X1() As Double, _
ByRef Y1() As Double, _
ByVal N As Long, _
ByVal DT As Double, _
ByRef PX() As Double, _
ByRef PY() As Double, _
ByRef PXY1() As Double, _
ByRef PXY2() As Double)
' Входные параметры:
' X1, Y1 - одномерные массивы длины N+1, первые N элементов которых
' содержат значения исходных реализаций случайных процессов
' Xa+s, Yb+s, s = 0, 1, ..., N-1 (тип: вещественный);
' N - количество заданных значений каждого из двух исходных рядов,
' N >= 4 - целая степень числа два (тип: целый);
' DT - заданный шаг сетки по оси времени, DT > 0 (тип: вещественный);
' Выходные параметры:
' PX - одномерный массив длины N/2+1, содержащий вычисленные значения
' автопериодограммы процесса X (тип: вещественный);
' PY - одномерный массив длины N/2+1, содержащий вычисленные значения
' автопериодограммы процесса Y (тип: вещественный);
' PXY1 - одномерный массив длины N/2+1, содержащий вычисленные значения
' вещественной части взаимной периодограммы процессов X, Y
' (тип: вещественный);
' PXY2 - одномерный массив длины N/2+1, содержащий вычисленные значения
' мнимой части взаимной периодограммы процессов X, Y (тип:
' вещественный).
Dim I As Long, JS As Long, NC As Long
Dim A1 As Double, A2 As Double, B1 As Double, B2 As Double
Dim G As Double, P As Double, Q As Double, R As Double
NC = N / 2 + 1
G = DT / (N * 4#)
P = 1#
Call FTF1C(X1, Y1, N, 1, 1, P)
X1(N + 1) = X1(1)
Y1(N + 1) = Y1(1)
For I = 1 To NC
JS = N + 2 - I
Q = X1(I)
R = X1(JS)
A1 = Q + R
B2 = R - Q
Q = Y1(I)
R = Y1(JS)
A2 = Q - R
B1 = Q + R
PXY1(I) = (A1 * B1 + A2 * B2) * G
PXY2(I) = (A1 * B2 - A2 * B1) * G
PY(I) = (B1 * B1 + B2 * B2) * G
PX(I) = (A1 * A1 + A2 * A2) * G
Next I
End Sub
' Вычисление дискретного или обратного дискретного преобразования Фурье
' комплексного ряда длины, равной степени двух, методом быстрого
' преобразования Фурье.
' Замечания по использованию:
' 1. Длина преобразуемого ряда L=N/K должна быть целой степенью двух:
' N/K = 2^M >= 4; 1 <= NE <= K.
' 2. Если преобразуемый ряд целиком заполняет массив A, то NE = K = 1.
' 3. Если вычислялось ДПФ, то на выходе P = 3.1415,
' а если OДПФ, то P = -3.1415.
' 4. Если на выходе P = 0, то элементы ряда X располагаются в двоично-
' инверсном порядке, а преобразование Фурье не выполняется.
' 5. Значения элементов массива A, не используемых при вычислении
' коэффициентов Фурье, при работе подпрограммы сохраняются.
Public Sub FTF1C(ByRef ARE() As Double, _
ByRef AIM() As Double, _
ByVal N As Long, _
ByVal NE As Long, _
ByVal K As Long, _
ByRef P As Double, _
Optional ByVal IFFT As Boolean = False)
' Параметры:
' ARE, AIM - вещественные одномерные массивы длины N,
' компоненты которых перед началом работы
' подпрограммы содержат сооответственно
' действительные и мнимые части элементов
' преобразуемого комплексного ряда X.
' По окончании работы подпрограммы содержат
' сооответственно действительные и мнимые части
' элементов преобразованного ряда Y на тех же
' местах, что и преобразуемый ряд X;
' P - заданная вещественная переменная, признак преобразования
' Фурье. При P > 0 выполняется дискретное преобразование
' Фурье, при P < 0 выполняется обратное дискретное
' преобразование Фурье.
' Входные параметры:
' N - заданная длина массива A (тип: целый);
' NE - заданный номеp первого элемента преобразуемого ряда
' в массиве A (тип: целый);
' K - заданный шаг выбора элементов преобразуемого ряда
' (тип: целый);
' IFFT - параметр нормировки значений обратного дискретного
' преобразования Фурье
Dim I As Long, II As Long, J As Long, JJ As Long
Dim N1 As Long, N2 As Long, M As Long, MM As Long
Dim R As Double, T As Double, SI As Double, CO As Double, S As Double
Dim C As Double, A As Double, B As Double, SYS001 As Double
SYS001 = 4 * Atn(1)
N1 = N / 2
MM = N1 / 2
N2 = N1 + K
J = NE
JJ = J
1 J = J + K
If J > N1 Then GoTo 7
II = JJ + N1
R = ARE(J)
ARE(J) = ARE(II)
ARE(II) = R
R = AIM(J)
AIM(J) = AIM(II)
AIM(II) = R
J = J + K
M = MM
3 If JJ <= M Then GoTo 5
JJ = JJ - M
M = M / 2
GoTo 3
5 JJ = JJ + M
If JJ <= J Then GoTo 1
R = ARE(J)
ARE(J) = ARE(JJ)
ARE(JJ) = R
R = AIM(J)
AIM(J) = AIM(JJ)
AIM(JJ) = R
I = J + N2
II = JJ + N2
R = ARE(I)
ARE(I) = ARE(II)
ARE(II) = R
R = AIM(I)
AIM(I) = AIM(II)
AIM(II) = R
GoTo 1
7 I = K
T = SYS001
If P > 0 Then
T = -T
ElseIf P = 0 Then
Exit Sub
End If
P = -T
10 SI = 0#
CO = 1#
S = Sin(T)
C = Cos(T)
T = 0.5 * T
II = I
I = I + I
For M = NE To II Step K
For J = M To N Step I
JJ = J + II
A = ARE(JJ)
B = AIM(JJ)
R = A * CO - B * SI
ARE(JJ) = ARE(J) - R
ARE(J) = ARE(J) + R
R = B * CO + A * SI
AIM(JJ) = AIM(J) - R
AIM(J) = AIM(J) + R
Next J
R = C * CO - S * SI
SI = C * SI + S * CO
CO = R
Next M
If I < N Then GoTo 10
If P < 0 And IFFT Then
For M = NE To N Step K
ARE(M) = ARE(M) / N * K
AIM(M) = AIM(M) / N * K
Next M
End If
End Sub