Привожу текст программы
| Код |
|---|
Function HILBERT(ByRef X() As Single, ByVal TS As Double, ByVal WMIN As Double, ByVal WMAX As Double) As Single()
Try
Dim DEN() As Complex = {Nothing, COMPLEXUNITY()}
Dim NUM(), Z() As Complex
Dim WP(2) As Double
Dim Y() As Single
Dim I, N As Integer
WP(1) = WMIN ' Нижняя частота среза
WP(2) = WMAX ' Верхняя частота среза
' Нормировать к частоте Найквиста
For I = 1 To 2
WP(I) = WP(I) * TS ' Конец полосы пропускания
Select Case WP(I)
Case Is <= 0
WP(I) = 0.01
Case Is > 0.499
WP(I) = 0.499
End Select
Next I
' Синтез нерекурсивного полосового фильтра с использованием весовых окон
NUM = Fir1(30, WP, Filtertype.BandPass)
' Фильтровать исходные данные
Y = filtfilt(NUM, DEN, X)
' Вычислить преобразование Гильберта
Z = Hilbert2(Y)
' Сформировать огибающую сигнала
N = Z.GetUpperBound(0)
ReDim Y(N)
For I = 1 To N
Y(I) = CSng(Rho(Z(I)))
Next I
Return Y
Catch E As Exception
Throw New ArgumentException(E.Message)
End Try
End Function
Function Hilbert2(ByVal X() As Single) As Complex()
Try
Dim J, K, M, N, NFFT As Integer
N = X.GetUpperBound(0)
Dim XSEG() As Complex = fft(X)
NFFT = XSEG.GetUpperBound(0)
Dim H(NFFT) As Integer
K = CInt(2 * Fix(N * 0.5))
If K = N Then
M = CInt(N * 0.5)
H(1) = 1
H(M + 1) = 1
Else
M = CInt((N + 1) * 0.5)
H(1) = 1
End If
For J = 2 To M
H(J) = 2
Next J
For J = 1 To NFFT
XSEG(J) *= H(J)
Next J
XSEG = fft(XSEG, IFFT:=True)
ReDim Preserve XSEG(N)
Return XSEG
Catch E As Exception
Throw New ArgumentException(E.Message)
End Try
End Function |
На рисунке () видно, что характер огибающей неверно отражает характер измеренных данных в диапазоне значений 90...130 Н и -90...-150 Н.
Вы не подскажите, как разработать алгоритм преобразования.