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.
The frequency response of a pendulum type seismic sensor is given
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
wp2 - wf2
The filter gain from input X to point Y2 is
Gb = G1*G2 / (1 - A1*G1 - A2*G1*G2)
(-1/w2) / (1 - jwf /(wQf) - (wf2/w2))
Let Df = 1 - jwf /(wQf) - (wf2/w2)
Then Gb = -1/(w2
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.
March 11, 2007