d10r6.frm

VB数值分析

VERSION 5.00
Begin VB.Form Form1 
   Caption         =   "Form1"
   ClientHeight    =   4920
   ClientLeft      =   60
   ClientTop       =   345
   ClientWidth     =   4755
   LinkTopic       =   "Form1"
   ScaleHeight     =   4920
   ScaleWidth      =   4755
   StartUpPosition =   3  'Windows Default
   Begin VB.CommandButton Command1 
      Caption         =   "Command1"
      Height          =   375
      Left            =   3240
      TabIndex        =   0
      Top             =   4200
      Width           =   1215
   End
End
Attribute VB_Name = "Form1"
Attribute VB_GlobalNameSpace = False
Attribute VB_Creatable = False
Attribute VB_PredeclaredId = True
Attribute VB_Exposed = False
Private Sub Command1_Click()
    'PROGRAM D10R6
    'Driver for routine CONVLV
    N = 16
    N2 = 32
    M = 9
    PI = 3.1415926
    Dim DATA(16), RESPNS(9), RESP(16), ANS(32)
    For I = 1 To N
        DATA(I) = 0#
        If I >= (N / 2 - N / 8) And I <= (N / 2 + N / 8) Then
            DATA(I) = 1#
        End If
    Next I
    Print
    For I = 1 To M
        RESPNS(I) = 0#
        If I > 2 And I < 7 Then RESPNS(I) = 1#
        RESP(I) = RESPNS(I)
    Next I
    ISIGN = 1
    Call CONVLV(DATA(), N, RESP(), M, ISIGN, ANS())
    'Compare with a direct convolution
    Print Tab(5); "I       CONVLV       Expected"
    For I = 1 To N
        CMP = 0#
        For J = 1 To M / 2
          CMP = CMP + DATA(((I - J - 1 + N) Mod N) + 1) * RESPNS(J + 1)
          CMP = CMP + DATA(((I + J - 1) Mod N) + 1) * RESPNS(M - J + 1)
        Next J
        CMP = CMP + DATA(I) * RESPNS(1)
        Print Tab(5); Format$(I, "###");
        Print Tab(12); Format$(ANS(I), "#.#####0");
        Print Tab(25); Format$(CMP, "#.#####0")
    Next I
End Sub
Sub CONVLV(DATA(), N, RESPNS(), M, ISIGN, ANS())
    Dim FFT(32)
    For I = 1 To CInt(M - 1) / 2
        RESPNS(N + 1 - I) = RESPNS(M + 1 - I)
    Next I
    For I = CInt(M + 3) / 2 To N - CInt(M - 1) / 2
        RESPNS(I) = 0#
    Next I
    Call TWOFFT(DATA(), RESPNS(), FFT(), ANS(), N)
    NO2 = CInt(N / 2)
    For I = 1 To NO2 + 1
        If ISIGN = 1 Then
          DUM = ANS(2 * I - 1)
          DUM1 = FFT(2 * I - 1) * DUM - FFT(2 * I) * ANS(2 * I)
          ANS(2 * I - 1) = DUM1 / NO2
          DUM2 = FFT(2 * I - 1) * ANS(2 * I) + FFT(2 * I) * DUM
          ANS(2 * I) = DUM2 / NO2
        ElseIf ISIGN = -1 Then
          If DUM = 0# And ANS(2 * I) = 0 Then
            Print "deconvolving at a response zero"
            Exit Sub
          End If
          ANS1 = FFT(2 * I - 1) * DUM + FFT(2 * I) * ANS(2 * I)
          DUM1 = DUM * DUM + ANS(2 * I) * ANS(2 * I)
          ANS(2 * I - 1) = ANS1 / DUM1 / NO2
          ANS1 = FFT(2 * I) * DUM - FFT(2 * I - 1) * ANS(2 * I)
          DUM2 = DUM * DUM + ANS(2 * I) * ANS(2 * I)
          ANS(2 * I) = ANS1 / DUM2 / NO2
        Else
          Print " no meaning for ISIGN"
        End If
    Next I
    ANS(2) = ANS(2 * NO2 + 1)
    Call REALFT(ANS(), NO2, -1)
    Erase FFT
End Sub
Sub TWOFFT(DATA1(), DATA2(), FFT1(), FFT2(), N)
    C1R = 0.5
    C1I = 0#
    C2R = 0#
    C2I = -0.5
    For J = 1 To N
        FFT1(2 * J - 1) = DATA1(J)
        FFT1(2 * J) = DATA2(J)
    Next J
    Call FOUR1(FFT1(), N, 1)
    FFT2(1) = FFT1(2)
    FFT2(2) = 0#
    FFT1(2) = 0#
    N2 = 2 * (N + 2)
    For J = 2 To N / 2 + 1
        J2 = 2 * J
        CONJR = FFT1(N2 - J2 - 1)
        CONJI = -FFT1(N2 - J2)
        H1R = C1R * (FFT1(J2 - 1) + CONJR) - C1I * (FFT1(J2) + CONJI)
        H1I = C1I * (FFT1(J2 - 1) + CONJR) + C1R * (FFT1(J2) + CONJI)
        H2R = C2R * (FFT1(J2 - 1) - CONJR) - C2I * (FFT1(J2) - CONJI)
        H2I = C2I * (FFT1(J2 - 1) - CONJR) + C2R * (FFT1(J2) - CONJI)
        FFT1(J2 - 1) = H1R
        FFT1(J2) = H1I
        FFT1(N2 - J2 - 1) = H1R
        FFT1(N2 - J2) = -H1I
        FFT2(J2 - 1) = H2R
        FFT2(J2) = H2I
        FFT2(N2 - J2 - 1) = H2R
        FFT2(N2 - J2) = -H2I
    Next J
End Sub
Sub REALFT(DATA(), N, ISIGN)
      THETA = 6.28318530717959 / 2# / N
      C1 = 0.5
      If ISIGN = 1 Then
          C2 = -0.5
          Call FOUR1(DATA(), N, 1)
      Else
          C2 = 0.5
          THETA = -THETA
      End If
      WPR = -2# * Sin(0.5 * THETA) ^ 2
      WPI = Sin(THETA)
      WR = 1# + WPR
      WI = WPI
      N2P3 = 2 * N + 3
      For I = 2 To N / 2 + 1
          I1 = 2 * I - 1
          I2 = I1 + 1
          I3 = N2P3 - I2
          I4 = I3 + 1
          WRS = CSng(WR)
          WIS = CSng(WI)
          H1R = C1 * (DATA(I1) + DATA(I3))
          H1I = C1 * (DATA(I2) - DATA(I4))
          H2R = -C2 * (DATA(I2) + DATA(I4))
          H2I = C2 * (DATA(I1) - DATA(I3))
          DATA(I1) = H1R + WRS * H2R - WIS * H2I
          DATA(I2) = H1I + WRS * H2I + WIS * H2R
          DATA(I3) = H1R - WRS * H2R + WIS * H2I
          DATA(I4) = -H1I + WRS * H2I + WIS * H2R
          WTEMP = WR
          WR = WR * WPR - WI * WRI + WR
          WI = WI * WPR + WTEMP * WPI + WI
      Next I
      If ISIGN = 1 Then
          H1R = DATA(1)
          DATA(1) = H1R + DATA(2)
          DATA(2) = H1R - DATA(2)
      Else
          H1R = DATA(1)
          DATA(1) = C1 * (H1R + DATA(2))
          DATA(2) = C1 * (H1R - DATA(2))
          Call FOUR1(DATA(), N, -1)
      End If
End Sub
Sub FOUR1(DATA(), NN, ISIGN)
      N = 2 * NN
      J = 1
      For I = 1 To N Step 2
          If J > I Then
              TEMPR = DATA(J)
              TEMPI = DATA(J + 1)
              DATA(J) = DATA(I)
              DATA(J + 1) = DATA(I + 1)
              DATA(I) = TEMPR
              DATA(I + 1) = TEMPI
          End If
          M = N / 2
          While M >= 2 And J > M
              J = J - M
              M = M / 2
          Wend
          J = J + M
      Next I
      MMAX = 2
      While N > MMAX
          ISTEP = 2 * MMAX
          THETA = 6.28318530717959 / (ISIGN * MMAX)
          WPR = -2# * Sin(0.5 * THETA) ^ 2
          WPI = Sin(THETA)
          WR = 1#
          WI = 0#
          For M = 1 To MMAX Step 2
              For I = M To N Step ISTEP
                  J = I + MMAX
                  TEMPR = CSng(WR) * DATA(J) - CSng(WI) * DATA(J + 1)
                  TEMPI = CSng(WR) * DATA(J + 1) + CSng(WI) * DATA(J)
                  DATA(J) = DATA(I) - TEMPR
                  DATA(J + 1) = DATA(I + 1) - TEMPI
                  DATA(I) = DATA(I) + TEMPR
                  DATA(I + 1) = DATA(I + 1) + TEMPI
              Next I
              WTEMP = WR
              WR = WR * WPR - WI * WPI + WR
              WI = WI * WPR + WTEMP * WPI + WI
          Next M
          MMAX = ISTEP
      Wend
End Sub