# This is a band pass filter based on a Rats program by Eduard Pelz,
# itself based on the paper "The Band Pass Filter" by Lawrence
# J. Christiano and Terry J. Fitzgerald (1999).
# Inputs: ($pu,$pl,@x) all in one array, where
# @x  - series of data 
# $pu - maximum period of oscillation of desired component (2<=pl<pu<infinity).
# $pl - minimum period of oscillation of desired component
# Output: @fx, array with filtered data
# Recommended: for quarterly data, $pl=6, $pu=32 returns component 
# with periods between 1.5 and 8 yrs, or for monthly data: $pl=2, $pu=24 
# returns component with all periods less than 2 yrs.
# Note: When feasible, we recommend dropping 2 years of data from the beginning
# and end of the filtered data series.  These data are relatively poorly estimated.

sub bpfilter{
    $pu=$_[0];
    $pl=$_[1];
    for($i=2;$i<=$#_;$i++){$x[$i-1]=$_[$i]}
#*******************************************************************
#*This section is designed to remove the drift from a time series using
#*the formula:  drift = [x(m) - x(startl)] / (m-1).
#*******************************************************************
    use constant PI => 4*atan2(1,1);
    $m=$#x;
    $undrift=1;
    $drift=0;
    if($undrift==1){$drift=($x[$m]-$x[1])/($m-1)}
    for($i=1;$i<=$m;$i++){$xun[$i]=$x[$i]-$drift*($i-1)}

#*******************************************************************
#*Create the ideal B's then construct the AA matrix.
#*******************************************************************

    $b = 2*PI/$pl;
    $a = 2*PI/$pu;
    $bnot = ($b-$a)/PI;
    $bhat = $bnot/2;
    for($i=1;$i<=$m;$i++){$bb[$i]=(sin($i*$b)-sin($i*$a))/($i*PI)}
    for($i=1;$i<$m;$i++){$bb2[$i+1]=$bb[$i]}
    $bb2[1]=$bnot;
    for($i=1;$i<=$m;$i++){
	$k=1;
	for($j=$i;$j<=$m;$j++){
	    $AA[$i][$j]=$bb2[$k];
	    $AA[$j][$i]=$bb2[$k];
	    $k++;
	}
    }
    $AA[1][1]=$bhat;
    $AA[$m][$m]=$bhat;
    for($i=1;$i<$m;$i++){
	$AA[$i+1][1]=$AA[$i][1]-$bb2[$i];
	$AA[$m-$i][$m]=$AA[$i][1]-$bb2[$i];
    }
    for($i=1;$i<=$m;$i++){
	$fx[$i]=0;
	for($j=1;$j<=$m;$j++){
	    $fx[$i]+=$AA[$i][$j]*$xun[$j];
	}
    }
    shift(@fx);
#    $fx[0]=sin(PI/2);
    return @fx;
}
1;
