'***************************************************************** '* 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