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.
fig:bills_explan


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

(1)
and
,
(2)
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

,
(3)
and
.
(4)
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
fig:one_sq_pulse



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.
fig:one_sq_pulse_freq_v1



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.
fig:one_sq_pulse_freq_v2


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

,
(5)
the analytic Fourier transform calculated with the definition matching that used by IDL's FFT is
.
(6)
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.
fig:one_sq_pul_freq_v2_ana


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

.
(7)
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
.
(8)
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.
fig:rc_in_out


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) fig:fft_sampling


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

,
(9)
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
.
(10)
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.
fig:two_sq_pulse_freq_v1



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.
fig:two_rou_pulse_freq_v1


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

.
(11)
The highest frequency component is set by the Nyquist critical frequency, given by
,
(12)
where D is the time period between successive samples. The frequency components between f1 and fc are just the integer multiples of f1:
(13)
(14)
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.
fig:two_sq_pulse_freq_v2



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.
fig:two_sq_pulse_freq_v3



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.
fig:two_rou_pulse_freq_v2


Compute Numerical Transfer Function

The simplest estimate for a transfer function by the division of two FFTs uses the expression

.
(18)
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
.
(19)
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.
fig:num_t_fun_freq_v1



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.
fig:num_t_fun_freq_v2


Compute Analytic Transfer Function

The analytic expression for the transfer function of an RC low-pass filter is

,
(20)
where the phase angle is
,
(21)
and the filter's critical frequency, fo, is a characteristic of the filter related to its time constant, t = RC, by
.
(22)
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.
fig:ana_t_fun_freq_v1


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

.
(23)
The expression for the output from an RC low pass filter given a three square pulse input is
.
(24)
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.
fig:rc_in_out_pred


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

,
(25)
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
.
(26)
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
.
(27)
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
(28)
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.
fig:predict_diff


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

,
(29)
,
(30)
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)

SEBAC Psychology Server

Send comments/criticisms/speculations to
Date Last Reviewed : August 19, 2002