From e0046bfb84d19271981757c633d3597d3c435691 Mon Sep 17 00:00:00 2001 From: Laurence Withers Date: Fri, 15 Apr 2011 18:56:19 +0000 Subject: [PATCH] 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. --- src/libiir/000_TopHeader.h | 1 + src/libiir/500_response.c | 71 +++++++++++++ src/libiir/500_response.h | 52 ++++++++++ src/tests/zplot_filter.c | 206 +++++++++++++++++++++++++++++++++++++ 4 files changed, 330 insertions(+) create mode 100644 src/libiir/500_response.c create mode 100644 src/libiir/500_response.h create mode 100644 src/tests/zplot_filter.c diff --git a/src/libiir/000_TopHeader.h b/src/libiir/000_TopHeader.h index b5a678e..1b04e1d 100644 --- a/src/libiir/000_TopHeader.h +++ b/src/libiir/000_TopHeader.h @@ -9,6 +9,7 @@ #define HEADER_libiir /* standard includes, or includes needed for type declarations */ +#include /* options for text editors kate: replace-trailing-space-save true; space-indent true; tab-width 4; diff --git a/src/libiir/500_response.c b/src/libiir/500_response.c new file mode 100644 index 0000000..b19ee9b --- /dev/null +++ b/src/libiir/500_response.c @@ -0,0 +1,71 @@ +/* libiir/src/libiir/500_response.c + * + * Copyright: ©2011, Laurence Withers. + * Author: Laurence Withers + * 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(jω) + * 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 +*/ diff --git a/src/libiir/500_response.h b/src/libiir/500_response.h new file mode 100644 index 0000000..5c4e62a --- /dev/null +++ b/src/libiir/500_response.h @@ -0,0 +1,52 @@ +/* libiir/src/libiir/500_response.h + * + * Copyright: ©2011, Laurence Withers. + * Author: Laurence Withers + * 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 +*/ diff --git a/src/tests/zplot_filter.c b/src/tests/zplot_filter.c new file mode 100644 index 0000000..d49ed9a --- /dev/null +++ b/src/tests/zplot_filter.c @@ -0,0 +1,206 @@ +/* libiir/src/tests/zplot_filter.c + * + * Copyright: ©2010–2011, Laurence Withers. + * Author: Laurence Withers + * License: GPLv3 +*/ + +#include "iir.h" + +#include +#include +#include +#include +#include +#include +#include +#include + + + +#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 +*/