Add ability to compute filter response

Adds functionality to compute the response of the filter's transfer function at
a given frequency. This uses the Z domain transform. A new test program
zplot_filter (which is compatible with the existing plot_filter) is provided.
This commit is contained in:
Laurence Withers 2011-04-15 18:56:19 +00:00
parent 4c0680fb02
commit e0046bfb84
4 changed files with 330 additions and 0 deletions

View File

@ -9,6 +9,7 @@
#define HEADER_libiir
/* standard includes, or includes needed for type declarations */
#include <complex.h>
/* options for text editors
kate: replace-trailing-space-save true; space-indent true; tab-width 4;

71
src/libiir/500_response.c Normal file
View File

@ -0,0 +1,71 @@
/* libiir/src/libiir/500_response.c
*
* Copyright: ©2011, Laurence Withers.
* Author: Laurence Withers <l@lwithers.me.uk>
* License: GPLv3
*/
static double complex
iir_response_aux(int nc, const double* c, int nd, const double* d, double freq)
{
int i;
double complex sum, resp, N;
/* Z transform; substitutions:
*
* z = exp()
* where j is the imaginary unit (written as I in the C code below)
* and where ω = 2πf
*
* z^-n = exp(jωn)
*
* exp(jx) = cos(x) + j.sin(x)
*
* z^-n = cos(-ωn) + j.sin(-ωn)
*
* N = -ω = -2πf
*/
N = -2 * M_PI * freq;
sum = 0;
for(i = 0; i < nc; ++i) {
sum += c[i] * (cos(i*N) + I * sin(i*N));
}
resp = sum;
sum = 1;
for(i = 1; i <= nd; ++i) {
sum += d[i - 1] * (cos(i*N) + I * sin(i*N));
}
return resp / sum;
}
double complex
iir_response_c(const struct iir_coeff_t* coeff, double freq)
{
return iir_response_aux(coeff->nc, coeff->c, coeff->nd, coeff->d, freq);
}
double complex
iir_response(const struct iir_filter_t* fi, double freq)
{
double complex resp;
for(resp = 1.0; fi; fi = fi->next) {
resp *= iir_response_aux(fi->nc, fi->c, fi->nd, fi->d, freq);
}
return resp;
}
/* options for text editors
vim: expandtab:ts=4:sw=4
*/

52
src/libiir/500_response.h Normal file
View File

@ -0,0 +1,52 @@
/* libiir/src/libiir/500_response.h
*
* Copyright: ©2011, Laurence Withers.
* Author: Laurence Withers <l@lwithers.me.uk>
* License: GPLv3
*/
/*! \defgroup response Measuring IIR response
These functions allow measuring the IIR response (i.e. frequency and phase of
the transfer function) at a given frequency. The response is a complex number;
see \c complex(5) for details. The magnitude \c cabs(z) is the gain of the
filter and the phase \c carg(z) its phase delay.
*/
/*!@{*/
/*! \brief Find response of coefficients at a single frequency.
\param coeff Set of coefficients.
\param freq Frequency of interest, as a fraction of sample rate.
\returns Response as complex number.
*/
double complex iir_response_c(const struct iir_coeff_t* coeff, double freq);
/*! \brief Find response of filter at a single frequency.
\param fi Filter to evaluate.
\param freq Frequency of interest, as a fraction of sample rate.
\returns Response as a complex number.
Similar to \ref iir_response_c(), except that it works on an instantiated
filter (and takes all parts of a chained filter into account). Does not alter
the state of the filter in any way and calls may be interleaved with filter
evaluation.
*/
double complex iir_response(const struct iir_filter_t* fi, double freq);
/*!@}*/
/* options for text editors
vim: expandtab:ts=4:sw=4:syntax=c.doxygen
*/

206
src/tests/zplot_filter.c Normal file
View File

@ -0,0 +1,206 @@
/* libiir/src/tests/zplot_filter.c
*
* Copyright: ©20102011, Laurence Withers.
* Author: Laurence Withers <l@lwithers.me.uk>
* License: GPLv3
*/
#include "iir.h"
#include <math.h>
#include <ctype.h>
#include <errno.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <unistd.h>
#include <complex.h>
#define NPOINTS (1000)
#define STEADY_STATE_CYCLES (20)
char* tmp_fname;
void
unlink_tmpfile(void)
{
unlink(tmp_fname);
}
int
do_plot(const char* filter_desc, double samp_rat, const char* png_filename)
{
int fd, ret;
FILE* fp;
char cmd_file[] = "/tmp/libiir-plot_filter.cmd.XXXXXX",
cmd[200];
fd = mkstemp(cmd_file);
if(fd == -1) {
perror("mkstemp");
return -1;
}
fp = fdopen(fd, "w");
fprintf(fp, "set terminal png size 1000,1000\n"
"set output '%s'\n"
"set multiplot layout 2,1 title \"Bode plot for filter '",
png_filename);
ret = 0;
while(*filter_desc) {
if(isspace(*filter_desc)) {
if(!ret) {
ret = 1;
putc('\n', fp);
}
} else {
putc(*filter_desc, fp);
}
++filter_desc;
}
fprintf(fp, "' at %fHz\"\n"
"set grid\n"
"set logscale\n"
"set ytics add ('-3dB' %f)\n"
"set xlabel 'Frequency (Hz)'\n"
"set ylabel 'Gain'\n"
"plot '%s' using 1:2 notitle\n"
"unset logscale y\n"
"set yrange [-180:180]\n"
"set ytics -180,45,180\n"
"set ylabel 'Phase (degrees)'\n"
"plot '%s' using 1:3 notitle\n",
samp_rat,
pow(10, -3.0/20),
tmp_fname,
tmp_fname);
fclose(fp);
snprintf(cmd, sizeof(cmd), "gnuplot %s", cmd_file);
ret = system(cmd);
unlink(cmd_file);
return ret;
}
double
interp(int step, int max, double start, double end)
{
return start + step * ((end - start) / (max - 1));
}
void
calc_response(FILE* fp,
const struct iir_filter_t* fi,
double samp_rat,
double start_freq,
double end_freq)
{
int step;
double freq;
double complex resp;
for(step = 0; step < NPOINTS; ++step) {
freq = exp(interp(step, NPOINTS, log(start_freq), log(end_freq)));
resp = iir_response(fi, freq / samp_rat);
fprintf(fp, "%e\t%e\t% 6.2f\n",
freq,
cabs(resp),
carg(resp) * 180 / M_PI);
}
}
int
safe_strtod(const char* str, double* d)
{
char* endp = 0;
errno = 0;
*d = strtod(str, &endp);
if(errno || !endp || *endp) return -1;
return 0;
}
int
main(int argc, char* argv[])
{
int fd;
double samp_rat, start_freq, end_freq;
FILE* fp;
struct iir_filter_t* fi;
/* process commandline arguments */
if(argc == 2 && !strcmp(argv[1], "--print-summary")) {
fputs("Generates Bode plot for a filter.\n", stdout);
return 0;
}
if(argc != 6) {
fputs("Usage: plot_filter 'filter_desc' samp_rat start_freq end_freq out.png\n",
stderr);
return 1;
}
fi = iir_parse(argv[1]);
if(!fi) {
fputs("Invalid filter description string.\n", stderr);
return 1;
}
if(safe_strtod(argv[2], &samp_rat) || samp_rat < 1e-6) {
fputs("Invalid sample rate. Positive float in Hz.\n", stderr);
return 1;
}
if(safe_strtod(argv[3], &start_freq) || start_freq < 1e-6) {
fputs("Invalid start frequency. Positive float in Hz.\n", stderr);
return 1;
}
if(safe_strtod(argv[4], &end_freq) || end_freq < 1e-6
|| end_freq > samp_rat || end_freq < start_freq)
{
fputs("Invalid end frequency. Positive float in Hz, less than sample\n"
"rate, but greater than start frequency.\n", stderr);
return 1;
}
printf("DC response: %f\n", cabs(iir_response(fi, 0)));
/* create temporary file for results; gnuplot will use this */
tmp_fname = strdup("/tmp/libiir-plot_filter.data.XXXXXX");
fd = mkstemp(tmp_fname);
if(fd == -1) {
perror("mkstemp");
return 1;
}
atexit(unlink_tmpfile);
fp = fdopen(fd, "w");
calc_response(fp, fi, samp_rat, start_freq, end_freq);
fclose(fp);
/* clean up (for valgrind) */
iir_filter_free(fi);
/* draw the plot */
return do_plot(argv[1], samp_rat, argv[5]);
}
/* options for text editors
kate: replace-trailing-space-save true; space-indent true; tab-width 4;
vim: expandtab:ts=4:sw=4
*/