'*****************************************************************
'* NUMERICAL DECONVOLUTION OF A SUSPENDED CAPTOR'S MASS ACCELE-  *
'* RATION RESPONSE SIGNAL  TO OBTAIN SPEED AT BASIS AND COMPUTE  *
'* THE ACCELERATION SHOCK SPECTRUM OF THIS SPEED WITH OR WITHOUT *
'* SIGNAL FILTERING OF HIGH FREQUENCIES                          *
'* ------------------------------------------------------------- *
'* The shock captor is considered as a damped elementary oscil-  *
'* lator of one degree of freedom (M,K) with a resonance fre-    *
'* quency F0 and a relative damping factor dzeta. Q is the       *
'* quality factor used to compute the shock spectrum.            *
'*                                                               *
'*                          Basic version by J-P Moreau, Paris   *
'*                                (in simple precision)          *
'*                                  (www.jpmoreau.fr)            *
'*                                                               *
'* SAMPLE RUN:                                                   *
'*                                                               *
'* Read data from file signal.dat                                *
'* Write results to file vsignal.lst.                            *
'*****************************************************************
DEFINT I-N

NMAX = 1024  'Maximum number of points of input signal

DIM ACC(NMAX), VIT(NMAX)
DIM S(800)

    CLS
    PRINT
    PRINT " DECONVOLUTION OF A SIGNAL F(T)"
    PRINT
    INPUT " Input data file name (without .dat): ", NOM$
    PRINT
    J = 0

    NOM2$ = "v" + NOM$ + ".lst"
    NOM1$ = NOM$ + ".dat"

    OPEN NOM2$ FOR OUTPUT AS #2
    OPEN NOM1$ FOR INPUT AS #1

' read title
    INPUT #1, TITLE$
    PRINT #2, "DECONVOLUTION OF SIGNAL:"
    PRINT #2, TITLE$
' print resume of program's functions
    PRINT #2,
    PRINT #2, " THIS PROGRAM ALLOWS TO:"
    PRINT #2,
    PRINT #2, " - READ MEASUREMENTS IN A DATA FILE,"
    PRINT #2, " - FILTER THESE MEASUREMENTS (OPTIONAL),"
    PRINT #2, " - RECONSTRUCT THE SPEED AT BASIS BY DECONVOLUTION,"
    PRINT #2, " - CALCULATE THE ACCELERATION SHOCK SPECTRUM"
    PRINT #2,

' Frequency of captor
    INPUT #1, F0
    PRINT #2, " FREQUENCY OF CAPTOR: "; F0; " HERTZ"
    PRINT #2,

' rod damping
    INPUT #1, DZETA
    PRINT #2, using " DAMPING OF CAPTOR: #.###"; DZETA
    PRINT #2,

' Value of Q for the shock spectrum
    INPUT #1, Q
    PRINT #2, " Q VALUE FOR THE SHOCK SPECTRUM: "; Q
    PRINT #2,

' Number of frequencies for the shock spectrum
    INPUT #1, NFREQ
    PRINT #2, " NUMBER OF FREQUENCIES FOR SPECTRUM: "; NFREQ
    PRINT #2,

' FREQUENCIES MINI & MAXI OF SPECTRUM
    INPUT #1, FINF
    INPUT #1, FSUP
    PRINT #2, " BEGIN FREQUENCY OF SPECTRUM: "; FINF
    PRINT #2, " END FREQUENCY OF SPECTRUM  : "; FSUP
    PRINT #2,

' READ MEASUREMENTS
    GOSUB 1000
    CLOSE #1
    FMAX = 0.5 / TS
    PRINT #2, using " NYQUIST FREQUENCY: ####.##"; FMAX

' LIMIT MAX FREQUENCY TO NYQUIST
    IF FSUP > FMAX THEN FSUP = FMAX

' CALL DECONVOLUTION
    GOSUB 2000

' print basis speed obtained by deconvolution
    PRINT #2,
    PRINT #2,
    PRINT #2, "   3RD COLUMN: SPEED OBTAINED BY DECONVOLUTION FROM"
    PRINT #2, "   UNFILTERED MASS ACCELERATION."
    PRINT #2, "   4TH COLUMN: SPEED OBTAINED BY DECONVOLUTION FROM"
    PRINT #2, "   FILTERED MASS ACCELERATION IF IT EXISTS (VOID HERE)."
    PRINT #2,
    PRINT #2,
    PRINT #2, "    TIME (S)    ACC M/S2     BASIS SPEED M/S"
    PRINT #2, "    --------   -----------   ---------------"
    format$ = "    ##.#####   ####.######     ####.######  "
    FOR I = 1 TO NDATA
      PRINT #2, USING format$; (I - 1) * TS; ACC(I); VIT(I)
    NEXT I

' Compute acceleration shock spectrum from speed
    GOSUB 3000


'Print acceleration shock spectrum of basis speed
    DF = (FSUP - FINF) / (NFREQ - 1)
    PRINT #2,
    PRINT #2,
    PRINT #2, " THE ACCELERATION SHOCK SPECTRUM  IS COMPUTED FROM"
    PRINT #2, " THE FILTERED ACCELERATION IF IT EXISTS, ELSE FROM"
    PRINT #2, " THE UNFILTERED ACCELERATION.'"
    PRINT #2,
    PRINT #2, " FREQUENCY STEP OF SHOCK SPECTRUM: "; DF
    PRINT #2,
    PRINT #2,
    PRINT #2, "   FREQUENCY HZ  SPECTRUM +     SPECTRUM -  IN M/S2"
    PRINT #2, "   ------------  -----------    -----------        "

    format$ = "     ####.##     ####.######    ####.######"
    FOR I = 1 TO NFREQ
      FR = FINF + (I - 1) * DF
      PRINT #2, USING format$; FR; S(2 * I - 1); S(2 * I)
    NEXT I

'Closing section
    PRINT #2,
    PRINT #2, " End of file "; NOM2$
    PRINT #2,
    CLOSE #2
    print " Program done."
	
END 'of main program


1000 'SUBROUTINE READATA(ACC,TS,NDATA)
'------------------------------------------------
'Read ndata: number of points of input signal
'        ts: sampling duration of signal in sec.
'    ACC(i): Table with ndata acceleration values
'------------------------------------------------ 
  INPUT #1, NDATA
  PRINT #2, NDATA; " POINTS READ."
  INPUT #1, TS
  PRINT #2, using " SAMPLING DURATION: #.######"; TS
  FOR I = 1 TO NDATA
    INPUT #1, TEMP, ACC(I)
  NEXT I
RETURN


2000 'SUBROUTINE DECON(ACC,TS,F0,DZETA,NDATA,VIT)
'--------------------------------------------------
' DECONVOLUTION SUBROUTINE
'
' Inputs:
'            ACC(I)  input mass acceleration
'            TS      sampling time of signal
'            F0      resonance frequency of captor
'            DZETA   relative damping factor
'            NDATA   number of points of ACC(I)
' Output:
'            VIT(I)  speed of basis by deconvolution
'---------------------------------------------------
  PI = 3.1415926535
  Q = 1 / (2*DZETA)
  OMEGA = 2 * PI * F0
  OMEGA2 = OMEGA * OMEGA
  OMEGAT = OMEGA * TS
  OMEGAQ = OMEGA * Q
  OMEGAQT = OMEGA * Q * TS
  EXPA = EXP(-OMEGAQT)
  T2SUR6 = TS * TS / 6
  Q0 = (1 - EXPA) / (OMEGA * OMEGAT) + TS / 2
  Q1 = TS + Q / OMEGA * EXPA - Q0
  Q2 = OMEGAQ * EXPA
  YN = 0
  XN = 0
  XPN = 0
  XPPN = ACC(1)
  VIT(1) = XPN
  FOR I = 2 TO NDATA
    YNM1 = YN
    XNM1 = XN
    XPNM1 = XPN
    XPPNM1 = XPPN
    XPPN = ACC(I)
    YPN = XPNM1 + Q0 * XPPN + Q1 * XPPNM1 + Q2 * (XNM1 - YNM1)
    XPN = XPNM1 + 0.5 * TS * (XPPN + XPPNM1)
    XN = XNM1 + TS * XPNM1 + T2SUR6 * (2 * XPPNM1 + XPPN)
    YN = XN + XPPN / OMEGA2 + (XPN - YPN) / OMEGAQ
    VIT(I) = YPN
  NEXT I
RETURN


3000 'SUBROUTINE SPEC(VIT,TS,FMIN,FMAX,DZETA,NDATA,TS,NFREQ,S)
'-------------------------------------------------------------
' ACCELERATION SHOCK SPECTRUM OF A SPEED
'
' Inputs:
'            VIT(I)  speed of basis from deconvolution
'            TS      sampling time of speed
'            FMIN    begin frequency of spectrum (<>0)
'            FMAX    end frequency of spectrum (<=Nyquist)
'            DZETA   relative damping factor
'            NDATA   number of points of ACC(I)
'            NFREQ   number of frequencies of spectrum (<=400)
' Output:
'            S(I)    shock spectrum (negative & positive)
'-------------------------------------------------------------
  DELTAF = (FSUP - FINF) / (NFREQ - 1)
  FMIN = FINF
' The algorithm fails if FMIN=0
  IF FMIN < DELTAF THEN FMIN = DELTAF
  FOR II = 1 TO NFREQ
    FR = FMIN + (II - 1) * DELTAF
    'CALL OSCIL(FR,VIT,ACC,T,DZETA,NDATA)
    GOSUB 4000
    'CALL EXTREM(ACC,YMIN,YMAX,NDATA)
    GOSUB 5000
    S(2 * II - 1) = YMAX
    S(2 * II) = YMIN
  NEXT II
RETURN


4000 'SUBROUTINE OSCIL(FREQU,VIT,ACC,TS,DZETA,NDATA)
'*************************************************************
'* This subroutine calculates the acceleration of the sismic *
'* mass of an elementary oscillator (1 degree of freedom),   *
'* the basis of which is submitted to a given speed, VIT(t). *
'* The input signal VIT is digitalized with a constant time  *
'* step, TS. The table VIT(I) contains NDATA speed values.   *
'* --------------------------------------------------------- *
'* INPUTS:                                                   *
'*         FREQU........: eigen frequency of oscillator      *
'*         VIT..........: input speed signal of basis VIT(t) *
'*         TS...........: time step of signal (constant)     *
'*         DZETA........: reduced damping factor             *
'*         NDATA........: number of points of signal         *
'* OUTPUT:                                                   *
'*         ACC..........: table containing acceleration res- *
'*                        ponse of oscillator (NDATA points) *
'*                                                           *
'*                      Basic version by J-P Moreau, Paris   *
'*************************************************************
  FREQU = FR: T = TS
  OMEGA = 2 * PI * FREQU
  OMEGA2 = OMEGA * OMEGA
  DELTA = SQR(1 - DZETA * DZETA)
  EXPA = EXP(-OMEGA * DZETA * T)
  GOMEGA = OMEGA * DELTA
  ARG = GOMEGA * T
  SINARG = SIN(ARG)
  COSARG = COS(ARG)
  QSI = DZETA / DELTA
  Q0 = (COSARG - QSI * SINARG) * EXPA
  Q1 = OMEGA2 * SINARG / GOMEGA * EXPA
  Q2 = (1# + (QSI * SINARG - COSARG) * EXPA) / T
  R1 = SINARG / GOMEGA * EXPA
  R0 = COSARG * EXPA + DZETA * OMEGA * R1
  R2 = (1 - DZETA * OMEGA * T) * R1 / T - COSARG * EXPA
  R3 = 1 - R1 / T
  XPNP1 = VIT(1)
  YPPNP1 = 0
  YPNP1 = 0
  ACC(1) = Q2 * XPNP1
  FOR I = 2 TO NDATA
    YPN = YPNP1
    YPPN = YPPNP1
    XPN = XPNP1
    XPNP1 = VIT(I)
    YPPNP1 = Q0 * YPPN + Q1 * (XPN - YPN) + Q2 * (XPNP1 - XPN)
    ACC(I) = YPPNP1
    YPNP1 = R0 * YPN + R1 * YPPN + R2 * XPN + R3 * XPNP1
  NEXT I
RETURN


5000 'SUBROUTINE EXTREM(ACC,YMIN,YMAX,NDATA)
'-----------------------------------------------------------
' SEEK MINIMUM AND MAXIMUM OF A TABLE ACC(I)
'
' Inputs:
'            ACC(I)  Table with NDATA values of acceleration
'            NDATA   Number of points of ACC(I)
' Outputs:
'            YMIN    minimum value of table ACC
'            YMAX    maximum value of table ACC
'-----------------------------------------------------------
  YMIN = 0: YMAX = 0
  FOR I = 1 TO NDATA
    IF ACC(I) > YMAX THEN YMAX = ACC(I)
    IF ACC(I) < YMIN THEN YMIN = ACC(I)
  NEXT I
RETURN

'End of file deconv.bas