Period-extending filter algorithm and analysis

March 10, 2007.

Please find attached a revised GIF sketch of my period extending digital filter and its pseudocode.

I have completed an analysis of this filter, and found it to be mathematically correct. My approach was to draw the analog equivalent of the filter and analyze its frequency response. I found that the sensor's response was exactly replaced by the filter's second order response. I have also tested the algorithm in an updated version of my SeisSim program, and it performed as expected. This filter compensates for under- or over- damping, and can make the output period longer (or shorter) than the sensor period.



Filter Analysis

  The frequency response of a pendulum type seismic sensor is given by
  Gp(w) = 1/ (1-(wp/ w)2) - j(wp/w)/Qp)
  The above equation holds for acceleration, velocity, and displacement, so long as we compare output with the same type of input.
  Let Dp = 1-(wp/w)2) - j(wp/w)/Qp
  Gp(w) = 1/Dp
  To change the apparent period of the sensor and possibly its damping, we would like a filter whose gain is
  Gf = Dp/Df
  where Df = 1-(wf/w)2) - j(wf/w)/Qf
  The overall response of the sensor would then be G0 = (1/Dp)*(Dp/Df) = 1/Df.
  If Qf = 0.707, then 1/Df has the same gain versus frequency as a 2nd order Butterworth highpass filter with corner frequency wf.
  Figure 2 is the analog schematic representation of the filter. It has two amplifiers G1 and G2, of gain = -j/w, which makes them actually integrators. There are two feedback paths of gain A1 and A2, respectively, and three forward output paths which are summed to constitute the output signal. The output components are (1), the original input signal, (2) B1 times Y1, and (3) B2 times Y2.
  The filter constants have the following values:
  A1 = -wf/Qf
  A2 = -wf2
  B1 = wp/Qp - wf/Qf
  B2 = wp2 - wf2
  The filter gain from input X to point Y2 is

  Gb = G1*G2 / (1 - A1*G1 - A2*G1*G2)
  Gb = (-1/w2) / (1 - jwf /(wQf) - (wf2/w2))
  Let Df = 1 - jwf /(wQf) - (wf2/w2)
  Then Gb = -1/(w2 Df).
  The gain at Y1 is
  Ga = Gb/G2 = (-1/(w2Df)) /(-j/w) = -j/(wDf)
  The gain at the output Y0 is
  Gf = 1 + B1*Ga + B2*Gb
  Gf = 1 - j(wp/Qp - wf /Qf)/(wDf) - (wp2 - wf2)/(w2Df)
After some manipulation, we obtain:

  Gf = (Df + Dp - Df)/Df = Dp/Df
which is the desired end result.

Bob McClure
March 11, 2007