RC Low Pass Filter as a
Test Case for
Behavioral Transfer Function Computations
Robert Kessel
Code 8123
Naval Research Laboratory
January 31, 2000
Abstract
An RC low pass filter's input and output voltages can serve as a useful test
case for behavioral transfer function data reduction processing. Analytic
expressions exist for all stages of the linear analysis and consequently, make
this filter particularly useful in debugging numerical software.
Key words: linear systems analysis, frequency domain, variable interval
schedules, step transition, transfer function, debugging numerical software.
RC Low Pass Filter as a
Test Case for
Behavioral Transfer Function Computations
When
debugging numerical software it is very nice to have a known analytic
test case. For the computations required by a linear analysis of steady-state
behavioral dynamics, the RC low pass filter can be used to provide a
particularly handy test case. The processing to first estimate a behavioral
transfer function from time series reinforcement and response rate data as
well as the use of transfer function to predict behavior is shown in
Figure 1. The test case simply replaces
the behavioral data appearing in both phases with appropriate input and
output pairs from an RC low pass filter.
Figure 1 A schematic
representation of the algorithm used to, first, extract a transfer function
and then, second, predict behavior to a new schedule. The upper left images
illustrate the behavioral output to a pattern of reinforcers. The middle left
depicts the conversion of that data to frequencies and the division of the
output by the input. The lower left image represents the transfer function.
The left side of the right half of the figure illustrates a schedule and its
frequency-domain representation. The lower portion of the right side of
the figure illustrates that data being multiplied by the transfer function to
produce the frequency representation of the predicted behavior. This is
transformed to the time domain and compared to the actual behavior as
illustrated in the upper right portion of the figure.
Because all the numerical routines used to generate the figures for this
discussion are written in
IDL,
the array index for an n element array
runs from 0 to n-1. Both FORTRAN and MATLAB have the array
index running from 1 to n. Please keep this in mind when porting
the code. One should also note that multiple definitions for the Fourier
transform exist that can alter the normalization or transform the results
to a complex conjugate.
Get the Fast Fourier Transform running
Before trying to do anything on the linear analysis, it is a good thing to
first show that the Fast Fourier Transform (FFT) routine works. A
standard FFT routine will convert an array of 2n
numbers into another array of 2n numbers where
n is some
integer.1
There are two important things to note. First, all examples in this
discussion will have 256 elements. Second, the time and frequency scales
have to be generated independently of the FFT routine.
In most general form, both input and output arrays in the pair are
complex. For the present case shown in
Figure 1, the time domain data is strictly real.
The output of an FFT given a strictly real input will be complex, but will
also have a symmetry property called Hermitian. A Hermitian complex
array has the real part of the negative frequency component equal to the
real part of the positive frequency component and the imaginary part of
the negative frequency component equal to the negative of imaginary part
of the positive frequency component, or as a pair of expressions
and
where h(f) is complex function of frequency. How very nice. What
one really needs is a test case for just the FFT routine before working on
the more general RC low pass filter test case.2
Since the test case inputs to follow use square pulses, it only makes sense
to use a square pulse to check the FFT.
Figure 2 shows a 256 element
array where the first 32 elements are 1.0 and the all other elements are
0.0. If the time interval for the 256 samples is 2,000 seconds,
this corresponds to the pulse's falling edge at 250 seconds. There are
some subtleties involved in how exactly one maps from real analog data to
this digital representation that will be addressed below. For the present, note
that plot's abscissa is just the array index number. Using this array as the input,
the real and imaginary parts shown in
Figure 3
are the complex output of an FFT routine. Also shown is the complex
conjugate which results when using the alternative definition of the
Fourier transforms employed Press, Flannery,
Teukolsky, & Vetterling (1992). (The alternative definition was
also used in the original FORTRAN-based data reduction code). Although both
the real and imaginary parts of the FFT have symmetries, the exact form seems
not the match the Hermitian symmetry given by
Equations 1 and
2. The apparent lack of the Hermitian
symmetry is caused by the wrapped form used by the IDL and almost all other
FFT routines (Press et al., 1992). How to
unwrap the frequency domain FFT output will be covered below in the context of
generating the time and frequency scales. Magnitude and phase plots for the FFT's complex output array are shown in
Figure 4. The conversion
of a complex FFT array element to its corresponding magnitude and phase is
done as
and
Note that when using Equation 4
in computer code, one has to include appropriate trapping logic for the
and
cases and set fh(f) to the correct
value as exceptions. Additionally, the computer code must return the correct
phase when either or both 
and 
are less than zero so
that full range
- p <
fh(f)
� p is used.
Figure 2 A 256 element square
pulse as an FFT test input. The values of the array are in
one_sq_pulse.dat
and the IDL routine used to generate both the values and the plot is
fig_one_sq_pulse.pro
Figure 3 The
real part (upper panel) and imaginary part (lower panel) of a 256 element
FFT output. Also shown as dashed line is the complex conjugate which
simply changes the sign of the imaginary part. The complex conjugate
form would be obtained with the Numerical Recipes FFT.
The values of the array are in
one_sq_pulse_freq.dat
and the IDL routine used to generate both the values and the plot is fig_one_sq_pulse_freq_v1.pro.
Figure 4 The
magnitude (upper panel) and phase (lower panel) of a 256 element FFT
output. The IDL routine used to generate the plot is
fig_one_sq_pulse_freq_v2.pro.
One final point to note before moving on to the RC low pass test case,
the FFT doesn't compute the continuous Fourier transform. For the most
part the differences are slight and generally confined to aliasing effects
at higher frequencies. Aliasing is a process in which the amplitudes from
higher frequency components are spuriously transferred down to
lower frequency components. See
Press et al. (1992)
for a longer discussion and the underlying mathematics. For the single
pulse given by
the analytic Fourier transform calculated with the definition matching
that used by IDL's FFT is
Figure 5 shows magnitude
and phase of Equation 6 evaluated at
each of the FFT frequencies in
Figure 4.3
While the analytic expression drops faster than the FFT, the curves have
essentially the same form.
Figure 5 The
magnitude (upper panel) and phase (lower panel) of an analytic Fourier
transform of a square pulse. The plot is normalized so the dc component's
real part is 1.0. The values of the array are in
one_sq_pul_freq_ana.dat
and the IDL routine used to generate both the values and the plot is
fig_one_sq_pul_freq_v2_ana.pro.
Measurement phase low pass input and output pair
By considering a pair of rectangular pulses as an input for the low-pass filter,
simple analytic expressions will exist for both input and output. The expression
for the input is
It is a good idea to have the pulse boundaries matching the time bins used to
sample Equation 7. Since there
will be 256 bins over the 2,000 seconds of the trial, these boundaries have a
spacing of D = 2000/256 or 7.81250 seconds. For
this measurement phase:
| t1 = | 218.750 |
| t2 = | 414.062 |
| t3 = | 656.250 |
| t4 = | 796.875
|
To place the pulse edges at some time bin boundaries it's occasionally
necessary to adjust the tns by
�0.001 to get the floating point logical equality
statements operate properly. The set above doesn't require this kind of tweaking
though. The output voltage generated by an RC low-pass filter in response to
the input square step given by
Equation 7 is
The input of Equation 7 and
resulting output of Equation 8
assuming the values for t1, t2, t3 and
t4 given above and that t = 35 are
shown in Figure 6.
Figure 6 The output of an RC
low pass filter (upper panel) that would be generated from a two pulse input
(lower panel). The values of the t, Vout(t), and
Vin(t) arrays are in
two_pulse.dat
and the IDL routines used to generate the values are
make_r_of_t.pro and
make_lo_pass_b_of_t_v4.pro. The IDL routine used to generate the plot is
plot_input_output.pro.
There are three analytic properties and one numerical aspect of
Equation 8 that must be kept in mind.
First, it is somewhat of an approximation since most ways of building and
measuring such a circuit would include a dc level in steady state. Provided the
pulses are in the first half of the trial and t is
kept small enough, Vout(t) will have returned effectively
to zero by the end of the 2,000 second trial and one can ignore the dc level.
Second, the voltage reached by the end of a pulse has to be carried as the starting
voltage for the decay in the next interval. Third, each pulse decays independently.
Hence, the expression gets a bit on the long side by the fifth interval. The evaluation
of Equation 8 at a set of 256 discrete
time samples is quite important. Figure 7
shows interrelationship between Vin(t),
Vout(t), and the time bins near the rising edge of the first pulse.
The input can be easily and nearly perfectly modeled by the flat discrete line
segments since the rising and falling pulse edges were selected to match the time
bin boundaries. Equation 8,
on the other hand, curves throughout the time bin. After trying both bin average
and snap shot sampling approaches, at least in this case, the bin average sampling
has better performance.
Figure 7 Vin(t)
and Vout(t) superposed on the time bins near the first pulse's
rising edge. The dashed lines show the bin sample average used for
Vout(t)
Convert to frequency domain with FFT
With the test case input and output pair in hand, it would be spiffy to use the FFT
to numerically convert the data shown in
Figure 6 to the frequency domain.
IDL's FFT assumes that the forward Fourier transform is
where h(f) is the frequency domain function equivalent to the time function
H(t). Press et al.'s (1992) definition for
the Fourier transform is complex conjugate of
Equation 9 given by
The Equation 10 definition is employed
in the FORTRAN version of the data reduction software
(Palya et al., 1996,
2000).
Figure 8 shows the real and
imaginary parts of vin(f) that result when the IDL FFT is
run on 256 sample array generated from
Equation 7.
Figure 9 shows the real
and imaginary parts of vout(f) that result when the IDL FFT
is run on 256 sample array generated from the output of an RC low-pass filter as
given by Equation 8. The RC
low-pass output is significantly greater than zero only for lower frequencies.
As expected, the higher frequency components are near-zero.
Figure 8 The real part
(upper panel) and imaginary part (lower panel) for the two square pulse test case
input converted to the frequency domain with an FFT. The values of the array
are in
two_sq_pulse_freq.dat and
the IDL routine used to generate both the values and the plot is
fig_two_sq_pulse_freq_v1.pro.
Figure 9 The real part
(upper panel) and imaginary part (lower panel) for an RC low pass filters output
from two square pulses converted to the frequency domain with an FFT. The values
of the array are in
two_rou_pulse_freq.dat
and the IDL routine used to generate both the values and the plot is
fig_two_rou_pulse_freq_v1.pro.
Plotting the Fourier Transforms
If you're accustomed to working with them, the wrapped frequency domain
form used by most FFT routines is natural enough. Still it's nice to unwrap
and plot the arrays as a function of frequency on occasion. We'll take a
small digression from the steps laid out in
Figure 1 to show how the unwrapping
is done. The FFT frequency scaling has to be generated first to make such
plots.
The frequency scale for an FFT is determined by the time period between
to successive samples and the total time interval for the sampling. In the
present case this total time interval is trial length. The lowest non-zero
frequency component is
The highest frequency component is set by the Nyquist critical frequency,
given by
where D is the time period between successive
samples. The frequency components between f1 and
fc are just the integer multiples of f1:
and so forth until fc is reached. An FFT also generates the
corresponding set of negative frequencies:
|
| |
| (15) |
| (16) |
. |
| (17) |
| |
|
For reasons concerning the internal workings of an FFT, the wrapped form
of the frequency domain array starts with the zero frequency component,
then all the positive components up to the Nyquist frequency component,
then continues from the most negative component, and finishes with the
f-1 as the last entry in the array.
Press et al. (1992) gives the
details. For a 256 element array using IDL's element index definition, the
structure is shown in
Table 1. Also shown in
Table 1 are the frequency values that
are obtained if one assumes that the trial length Ttrial is
2,000 second and the sampling interval D
is 7.8125 seconds.
Table 1 Array indices and
corresponding frequencies output by the IDL FFT in wrapped form
| index | 0 | 1 |
2 | ... | 126
| 127 | 128 |
129 | ... | 255 |
| frequency | 0 |
f1 | f2 | ... |
(fc - 2 f1) |
(fc - f1) |
fc | -(fc - f1) | ...
| f-1 |
| (in mHz) | 0.0 | 0.5 |
1.0 | ... | 63.0
| 63.4 | 64.0 |
-63.5 | ... |
-0.5 |
Unwrapping the structure of Table 1 so
the FFT is in the order of increasing frequency involves swapping the two
halves of the array. Figure 10
shows the unwrapped two square pulse input as real and imaginary part plots.
The Hermitian symmetries given by
Equations 1
and 2 can be clearly seen.
Figure 11 and
Figure 12 show magnitude and
phase plots for the two pulse input and the RC low pass filter output.
Figure 10 The real part
(upper panel) and imaginary part (lower panel) for the two square pulse input
converted to the frequency domain with an FFT and then unwrapped in
frequency. The values of the f and vin(f) arrays are in
two_sq_pulse_freq_v2.dat
and the IDL routine used to generate both the values and the plot is
fig_two_sq_pulse_freq_v2.pro.
Figure 11 The
magnitude (upper panel) and phase (lower panel) for the two square pulse
input converted to the frequency domain with an FFT and then unwrapped in
frequency. The IDL routine used to generate the plot is
fig_two_sq_pulse_freq_v3.pro.
Figure 12 The
magnitude (upper panel) and phase (lower panel) for an RC low pass filters
output from two square pulses converted to the frequency domain with an
FFT. The values of the f and vout(f) arrays are in
two_rou_pulse_freq_v2.dat and the IDL routine used to generate
both the values and the plot is
fig_two_rou_pulse_freq_v2.pro.
Compute Numerical Transfer Function
The simplest estimate for a transfer function by the division of two FFTs uses
the expression
For the RC low pass test case, the output voltage, vout(f),
takes the place of b(f) and similarly the input voltage,
vin(f), replaces r(f) in
Equation 18 yielding a numerical estimate for
the RC low pass transfer function given by
Since vout(f) and vin(f) are complex,
complex arithmetic is needed within the computer program evaluating
Equation 19.
Figure 13 and
Figure 14
show the transfer function one obtains with
Equation 19 as real and imaginary
part or magnitude and phase plot pairs in unwrapped form.
Figure 13 The real part
(upper panel) and imaginary part (lower panel) for the transfer function
unwrapped in frequency. The FFT wrapped array of the transfer function alone
is in
num_t_fun_freq_v1.dat.
The values of the f and gRC,num(f) in unwrapped
arrays are in
num_t_fun_freq_v2.dat
and the IDL routine used to generate both the values and the plot is
fig_num_t_fun_freq_v1.pro.
Figure 14 The magnitude
(upper panel) and phase (lower panel) for the transfer function unwrapped in
frequency. The IDL routine used to generate the plot is
fig_num_t_fun_freq_v2.pro.
Compute Analytic Transfer Function
The analytic expression for the transfer function of an RC low-pass filter is
where the phase angle is
and the filter's critical frequency, fo, is a characteristic of
the filter related to its time constant, t = RC, by
Equation 20 is the complex
conjugate of the expression in
Palya et al., (1996)
corresponding to the switch in Fourier transform definition.
Figure 15 is the magnitude
and phase plots from direct evaluation of gRC(f).
Figure 14 agrees closely
with Figure 15.
Figure 15 The magnitude
(upper panel) and phase (lower panel) for the transfer function unwrapped in
frequency. The FFT wrapped array of the transfer function alone is in
ana_t_fun_freq_v1.dat.
The values of the f and gRC(f) in unwrapped
arrays are in
ana_t_fun_freq_v2.dat and
the IDL routine used to generate both the values and the plot is
fig_ana_t_fun_freq_v2.pro.
A rather interesting difficulty arose while developing the ILD routines
that have the expected close agreement between
Figure 14
and Figure 15. In the
initial versions of the two pulse input/output test pair, the rising and falling
edges had different time relationships. The output effectively lagged
the input for the rising edge, but led the input for the falling edge. This
kind of process is not describable within a linear analysis which
assumes the same time relationship on both edges. This aspect of the
processing probably bares further study, since the numerical estimate
for the RC low pass transfer function with this time relationship
mismatch looked more than a bit reminiscent of those seen with
behavioral data. Perhaps the birds employ two different processes
when adapting to rising and falling reinforcement rates.
Test phase low pass input and output pair
As the test phase of the linear analysis prediction consider a three
pulse input. The analytic form, analogous to
Equation 7, for the input is
The expression for the output from an RC low pass filter given a three
square pulse input is
For the test phase it is again useful to select the pulse edges that match
time bin boundaries:
| t1 = | 125.000 |
| t2 = | 195.312 |
| t3 = | 289.062 |
| t4 = | 359.375 |
| t5 = | 539.062 |
| t6 = | 695.312
|
The input of Equation 23
and resulting output of
Equation 24 assuming the
values for t1, t2, t3, t4,
t5 and t6 given above and that
t = 35 are shown in
Figure 16.
Figure 16 The output of an RC
low pass filter Vout(t) (upper panel) that would be generated from
a three pulse input Vin(t) (lower panel). Also plotted in the top panel,
though all but impossible to distinguish at this plot's scale, is the linear predict
Voutpred(t). The values of
the t, Vout(t), Vin(t), and
Voutpred(t) arrays are in
three_pulse.dat. The IDL routines
used to generate the analytic test phase values are
make_r_of_t.pro and
make_lo_pass_b_of_t_v5.pro. The linear predict is done with
IDL routine linear_predict.pro.
The IDL routine used to generate the plot is
plot_input_output_v2.pro.
Compute linear prediction
The prediction of the RC low pass filter's output based on the numerical
estimate of the transfer function in
Figure 1 given a three pulse input is
straightforward. One just follows the steps
shown on the right side of Figure 1.
The frequency domain plots obtained as intermediate results of the prediction
are very similar to those obtained while estimating the transfer function.
Consequently, these plots are not included.
The prediction requires three steps. First, the input square pulse
data is converted into the frequency domain with an FFT. Then one
uses the expression for the linear prediction in frequency domain
given by
where b(f), g(f), and r(f) are the Fourier transforms of
B(t), G(t-t�), and R(t�)
respectively. Substituting in the appropriate functions for the RC low pass
filter the expression becomes
As with the estimation of the transfer function, evaluation of
Equation 26 is done with complex
arithmetic. This leaves only the back transform of
voutpred(f) to the time domain. The IDL
definition of the back or inverse FFT is
Note that the IDL FFT routine also automatically handles the normalization of
the back transform. In contrast,
(Press et al., 1992) define the back
transform as
and the user is responsible for setting the normalization properly. The prediction
resulting from these three steps is shown in
Figure 16.
Compare results qualitatively
The qualitative comparison between the linear prediction and the analytic
expression is just a matter of looking at the two curves. As
noted in Figure 16, the linear
prediction is all but indistinguishable from the analytic expression given by
Equation 24. The residuals
plot in Figure 17 shows that the
prediction is accurate to the machine precision scale.
Figure 17 Residual for the
linear prediction Voutpred(t) - Vout(t)
as a function of time. The IDL routine used to generate the plot is
plot_output_diff.pro.
Compare results quantitatively
Absent an established reason to do otherwise, the reduced chi-squared sum of
squared discrepancies is the standard method to answer this question
(Bevington, 1969). Following
Bevington (1969), as well as
Palya et al. (1996),
cn2 for the discrepancy between measured
response rates B(t) and the predicted response rates
Bpred(t) is
where n is the difference between the number of data
points N and the number of adjustable parameters n and
s2Bi is
the variance in the ith measured response rate
B(ti). Since the predictions have zero free parameters, n = 0,
the normalization factor simplifies, as shown, to 1/N and
n = N.
Because the RC low pass test case is entirely a numerical exercise, there is not a
set of s2
Voutpred,i values
to use with either Equation 29 or
30. One is, of course, free to compute an
equal weighted
cn
2. Beyond establishing quantitatively that the prediction is
at the machine precision scale, little will be gained from such an exercise. A set of
s2Voutpred,i
values could be generated independently and then varied to gain a sense of how
the probabilities depends on the linear prediction's fidelity to the analytic
expression.
References
Bevington, P.R. (1969). Data reduction
and error analysis for the physical sciences (pp. 187-203). New York:
McGraw-Hill Book Company.
Palya, W.L., Walter, D., Kessel, R., and Lucke, R.
(1996). Investigating Dynamic Behavior with a Fixed-Time Extinction Schedule
and Linear Analysis. The Journal of the Experimental Analysis of
Behavior, 66, 391--409.
Palya, W.L., Walter, D., Kessel, R., and Lucke,
R. (1999). Linear Modeling of Steady-State Behavioral Dynamics. In press at
The Journal of the Experimental Analysis of Behavior.
Press, W.H., Flannery, B.P., Teukolsky, S.A., and
Vetterling, W.T. (1992). Numerical Recipes in FORTRAN 77: The Art of
Scientific Computing, 2nd (pp. 490-602). Cambridge:
Cambridge University Press.
Footnotes:
1 More specialized FFT packages do exist that handle
arrays of any length.
2 In
other words, a pre-test case before the test case. Fortunately, the sequence truncates
here. A pre-pre-test case before the pre-test case isn't required.
3 Ok,
the use of Equation 6 really does seem to
be a pre-pre-test case, but the sequence really does stop here. Really. ;-)
File translated from
LATEX
by TTH, version 2.21.
On 24 May 2000, 15:26.
(with later hand modifications for the equations, figure cross
references, and such)
Date Last Reviewed : August 19, 2002
,
,
.
,
.
.
.
,
.
.
,
.
.
,
,
.
.
.
,
.
.
, 
,