/*
* filterkit.c (library "filterkit.a"): Kaiser-windowed low-pass filter support.
*
/* filterkit.c
*
* LpFilter() - Calculates the filter coeffs for a Kaiser-windowed low-pass
* filter with a given roll-off frequency. These coeffs
* are stored into a array of doubles.
* writeFilter() - Writes a filter to a file.
* makeFilter() - Calls LpFilter() to create a filter, then scales the double
* coeffs into an array of half words.
* readFilter() - Reads a filter from a file.
* FilterUp() - Applies a filter to a given sample when up-converting.
* FilterUD() - Applies a filter to a given sample when up- or down-
* converting.
* initZerox() - Initialization routine for the zerox() function. Must
* be called before zerox() is called. This routine loads
* the correct filter so zerox() can use it.
* zerox() - Given a pointer into a sample, finds a zero-crossing on the
* interval [pointer-1:pointer+2] by iteration.
* Query() - Ask the user for a yes/no question with prompt, default,
* and optional help.
* GetUShort() - Ask the user for a unsigned short with prompt, default,
* and optional help.
* GetDouble() - Ask the user for a double with prompt, default, and
* optional help.
* GetString() - Ask the user for a string with prompt, default, and
* optional help.
*/
#import "resample.h"
#import "filterkit.h"
#import <libc.h>
#import <math.h>
/*
* Getstr() will print the passed "prompt" as a message to the user, and
* wait for the user to type an input string. The string is
* then copied into "answer". If the user types just a carriage
* return, then the string "defalt" is copied into "answer".
* ???
* "Answer" and "defalt" may be the same string, in which case,
* the default value will be the original contents of "answer".
* ???
*/
static void getstr(char *prompt, char *defalt, char *answer)
{
char *p,s[200];
printf("%s (default = %s):",prompt,defalt);
strcpy(s, prompt);
if (!(*s))
strcpy(s, "input:");
p = s;
while (*p && *p != '(')
p++;
p--;
while (*p == ' ')
p--;
*(++p) = '\0';
gets(answer);
if (answer[0] == '\0') {
strcpy(answer, defalt);
printf("\t%s = %s\n",s,answer);
} else {
printf("\t%s set to %s\n",s,answer);
}
}
/* LpFilter()
*
* reference: "Digital Filters, 2nd edition"
* R.W. Hamming, pp. 178-179
*
* Izero() computes the 0th order modified bessel function of the first kind.
* (Needed to compute Kaiser window).
*
* LpFilter() computes the coeffs of a Kaiser-windowed low pass filter with
* the following characteristics:
*
* c[] = array in which to store computed coeffs
* frq = roll-off frequency of filter
* N = Half the window length in number of coeffs
* Beta = parameter of Kaiser window
* Num = number of coeffs before 1/frq
*
* Beta trades the rejection of the lowpass filter against the transition
* width from passband to stopband. Larger Beta means a slower
* transition and greater stopband rejection. See Rabiner and Gold
* (Theory and Application of DSP) under Kaiser windows for more about
* Beta. The following table from Rabiner and Gold gives some feel
* for the effect of Beta:
*
* All ripples in dB, width of transition band = D*N where N = window length
*
* BETA D PB RIP SB RIP
* 2.120 1.50 +-0.27 -30
* 3.384 2.23 0.0864 -40
* 4.538 2.93 0.0274 -50
* 5.658 3.62 0.00868 -60
* 6.764 4.32 0.00275 -70
* 7.865 5.0 0.000868 -80
* 8.960 5.7 0.000275 -90
* 10.056 6.4 0.000087 -100
*/
#define IzeroEPSILON 1E-21 /* Max error acceptable in Izero */
static double Izero(double x)
{
double sum, u, halfx, temp;
int n;
sum = u = n = 1;
halfx = x/2.0;
do {
temp = halfx/(double)n;
n += 1;
temp *= temp;
u *= temp;
sum += u;
} while (u >= IzeroEPSILON*sum);
return(sum);
}
void LpFilter(double c[], int N, double frq, double Beta, int Num)
{
double IBeta, temp, inm1;
int i;
/* Calculate ideal lowpass filter impulse response coefficients: */
c[0] = 2.0*frq;
for (i=1; i<N; i++) {
temp = PI*(double)i/(double)Num;
c[i] = sin(2.0*temp*frq)/temp; /* Analog sinc function, cutoff = frq */
}
/*
* Calculate and Apply Kaiser window to ideal lowpass filter.
* Note: last window value is IBeta which is NOT zero.
* You're supposed to really truncate the window here, not ramp
* it to zero. This helps reduce the first sidelobe.
*/
IBeta = 1.0/Izero(Beta);
inm1 = 1.0/((double)(N-1));
for (i=1; i<N; i++) {
temp = (double)i * inm1;
c[i] *= Izero(Beta*sqrt(1.0-temp*temp)) * IBeta;
}
}
/* Write a filter to a file
* Filter file format:
* file name: "F" Nmult "T" Nhc ".filter"
* 1st line: the string "ScaleFactor" followed by its value.
* 2nd line: the string "Length" followed by Nwing's value.
* 3rd line: the string "Nmult" followed by Nmult's value.
* 4th line: the string "Coeffs:" on a separate line.
* following lines: Nwing number of 16-bit impulse response values
* in the right wing of the impulse response (the Imp[] array).
* (Nwing is equal to Npc*(Nmult+1)/2+1, where Npc is defined in the
* file "resample.h".) Each coefficient is on a separate line.
* next line: the string "Differences:" on a separate line.
* following lines: Nwing number of 16-bit impulse-response
* successive differences: ImpD[i] = Imp[i+1] - Imp[i].
* ERROR codes:
* 0 - no error
* 1 - could not open file
*/
int writeFilter(HWORD Imp[], HWORD ImpD[], UHWORD LpScl, UHWORD Nmult,
UHWORD Nwing)
{
char fname[30];
FILE *fp;
int i;
sprintf(fname, "F%dT%d.filter", Nmult, Nhc);
fp = fopen(fname, "w");
if (!fp)
return(1);
fprintf(fp, "ScaleFactor %d\n", LpScl);
fprintf(fp, "Length %d\n", Nwing);
fprintf(fp, "Nmult %d\n", Nmult);
fprintf(fp, "Coeffs:\n");
for (i=0; i<Nwing; i++) /* Put array of 16-bit filter coefficients */
fprintf(fp, "%d\n", Imp[i]);
fprintf(fp, "Differences:\n");
for (i=0; i<Nwing; i++) /* Put array of 16-bit filter coeff differences */
fprintf(fp, "%d\n", ImpD[i]);
fclose(fp);
printf("Wrote filter file '%s' in current directory.\n",fname);
return(0);
}
/* ERROR return codes:
* 0 - no error
* 1 - Nwing too large (Nwing is > MAXNWING)
* 2 - Froll is not in interval [0:1)
* 3 - Beta is < 1.0
* 4 - LpScl will not fit in 16-bits
*
* Made following global to avoid stack problems in Sun3 compilation: */
#define MAXNWING 8192
static double ImpR[MAXNWING];
int makeFilter(HWORD Imp[], HWORD ImpD[], UHWORD *LpScl, UHWORD Nwing,
double Froll, double Beta)
{
double DCgain, Scl, Maxh;
HWORD Dh;
int i, temp;
if (Nwing > MAXNWING) /* Check for valid parameters */
return(1);
if ((Froll<=0) || (Froll>1))
return(2);
if (Beta < 1)
return(3);
/*
* Design Kaiser-windowed sinc-function low-pass filter
*/
LpFilter(ImpR, (int)Nwing, 0.5*Froll, Beta, Npc);
/* Compute the DC gain of the lowpass filter, and its maximum coefficient
* magnitude. Scale the coefficients so that the maximum coeffiecient just
* fits in Nh-bit fixed-point, and compute LpScl as the NLpScl-bit (signed)
* scale factor which when multiplied by the output of the lowpass filter
* gives unity gain. */
DCgain = 0;
Dh = Npc; /* Filter sampling period for