libctsim/filter.cpp

   1 /*****************************************************************************
   2 ** FILE IDENTIFICATION
   3 **
   4 **     Name:                   filter.cpp
   5 **     Purpose:                Routines for signal-procesing filters
   6 **     Progammer:              Kevin Rosenberg
   7 **     Date Started:           Aug 1984
   8 **
   9 **  This is part of the CTSim program
  10 **  Copyright (C) 1983-2000 Kevin Rosenberg
  11 **
  12 **  $Id: filter.cpp,v 1.1 2000/06/19 02:59:34 kevin Exp $
  13 **
  14 **  This program is free software; you can redistribute it and/or modify
  15 **  it under the terms of the GNU General Public License (version 2) as
  16 **  published by the Free Software Foundation.
  17 **
  18 **  This program is distributed in the hope that it will be useful,
  19 **  but WITHOUT ANY WARRANTY; without even the implied warranty of
  20 **  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  21 **  GNU General Public License for more details.
  22 **
  23 **  You should have received a copy of the GNU General Public License
  24 **  along with this program; if not, write to the Free Software
  25 **  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
  26 ******************************************************************************/
  27
  28 #include "ct.h"
  29
  30
  31 /* NAME
  32  *   filter_generate                            Generate a filter
  33  *
  34  * SYNOPSIS
  35  *   f = filter_generate (filt_type, bw, xmin, xmax, n, param, domain, analytic)
  36  *   double f                           Generated filter vector
  37  *   int filt_type                      Type of filter wanted
  38  *   double bw                          Bandwidth of filter
  39  *   double xmin, xmax                  Filter limits
  40  *   int n                              Number of points in filter
  41  *   double param                       General input parameter to filters
  42  *   int domain                         FREQ or SPATIAL domain wanted
  43  *   int numint                         Number if intervals for calculating
  44  *                                      discrete inverse fourier xform
  45  *                                      for spatial domain filters.  For
  46  *                                      ANALYTIC solutions, use numint = 0
  47  */
  48
  49 double *
  50 filter_generate (const FilterType filt_type, double bw, double xmin, double xmax, int n, double param, const DomainType domain, int numint)
  51 {
  52   double *f = new double [n];
  53   double xinc = (xmax - xmin) / (n - 1);
  54
  55   if (filt_type == FILTER_SHEPP) {
  56     double a = 2 * bw;
  57     double c = - 4. / (a * a);
  58     int center = (n - 1) / 2;
  59     int sidelen = center;
  60     f[center] = 4. / (a * a);
  61
  62     for (int i = 1; i <= sidelen; i++ )
  63       f [center + i] = f [center - i] = c / (4 * (i * i) - 1);
  64   } else if (domain == D_FREQ) {
  65     double x;
  66     int i;
  67     for (x = xmin, i = 0; i < n; x += xinc, i++)
  68       f[i] = filter_frequency_response (filt_type, x, bw, param);
  69   } else if (domain == D_SPATIAL) {
  70     double x;
  71     int i;
  72     for (x = xmin, i = 0; i < n; x += xinc, i++)
  73       if (numint == 0)
  74         f[i] = filter_spatial_response_analytic (filt_type, x, bw, param);
  75       else
  76         f[i] = filter_spatial_response_calc (filt_type, x, bw, param, numint);
  77   } else {
  78     sys_error (ERR_WARNING, "Illegal domain %d [filt_generate]", domain);
  79     return (NULL);
  80   }
  81
  82   return (f);
  83 }
  84
  85
  86 /* NAME
  87  *   filter_spatial_response_calc       Calculate filter by discrete inverse fourier
  88  *                                      transform of filters's frequency
  89  *                                      response
  90  *
  91  * SYNOPSIS
  92  *   y = filter_spatial_response_calc (filt_type, x, bw, param, n)
  93  *   double y                   Filter's response in spatial domain
  94  *   int filt_type              Type of filter (definitions in ct.h)
  95  *   double x                   Spatial position to evaluate filter
  96  *   double bw                  Bandwidth of window
  97  *   double param               General parameter for various filters
  98  *   int n                      Number of points to calculate integrations
  99  */
 100
 101 double
 102 filter_spatial_response_calc (int filt_type, double x, double bw, double param, int n)
 103 {
 104   double zmin, zmax;
 105
 106   if (filt_type == FILTER_TRIANGLE) {
 107     zmin = 0;
 108     zmax = bw;
 109   } else {
 110     zmin = 0;
 111     zmax = bw / 2;
 112   }
 113   double zinc = (zmax - zmin) / (n - 1);
 114
 115   double z = zmin;
 116   double q [n];
 117   for (int i = 0; i < n; i++, z += zinc)
 118     q[i] = filter_frequency_response (filt_type, z, bw, param) * cos (TWOPI * z * x);
 119
 120   double y = 2 * integrateSimpson (zmin, zmax, q, n);
 121
 122   return (y);
 123 }
 124
 125
 126 /* NAME
 127  *    filter_frequency_response                 Return filter frequency response
 128  *
 129  * SYNOPSIS
 130  *    h = filter_frequency_response (filt_type, u, bw, param)
 131  *    double h                  Filters frequency response at u
 132  *    int filt_type             Type of filter
 133  *    double u                  Frequency to evaluate filter at
 134  *    double bw                 Bandwidth of filter
 135  *    double param              General input parameter for various filters
 136  */
 137
 138 double
 139 filter_frequency_response (int filt_type, double u, double bw, double param)
 140 {
 141   double q;
 142   double au = fabs (u);
 143
 144   switch (filt_type) {
 145   case FILTER_BANDLIMIT:
 146     if (au >= bw / 2)
 147       q = 0.;
 148     else
 149       q = 1;
 150     break;
 151   case FILTER_ABS_BANDLIMIT:
 152     if (au >= bw / 2)
 153       q = 0.;
 154     else
 155       q = au;
 156     break;
 157   case FILTER_TRIANGLE:
 158     if (au >= bw)
 159       q = 0;
 160     else
 161       q = 1 - au / bw;
 162     break;
 163   case FILTER_COSINE:
 164     if (au >= bw / 2)
 165       q = 0;
 166     else
 167       q = cos(PI * u / bw);
 168     break;
 169   case FILTER_ABS_COSINE:
 170     if (au >= bw / 2)
 171       q = 0;
 172     else
 173       q = au * cos(PI * u / bw);
 174     break;
 175   case FILTER_SINC:
 176     q = bw * sinc (PI * bw * u, 1.);
 177     break;
 178   case FILTER_ABS_SINC:
 179     q = au * bw * sinc (PI * bw * u, 1.);
 180     break;
 181   case FILTER_G_HAMMING:
 182     if (au >= bw / 2)
 183       q = 0;
 184     else
 185       q = param + (1 - param) * cos (TWOPI * u / bw);
 186     break;
 187   case FILTER_ABS_G_HAMMING:
 188     if (au >= bw / 2)
 189       q = 0;
 190     else
 191       q = au * (param + (1 - param) * cos(TWOPI * u / bw));
 192     break;
 193   default:
 194     q = 0;
 195     sys_error (ERR_WARNING,
 196                "Frequency response for filter %d not implemented [filter_frequency_response]",
 197                filt_type);
 198     break;
 199   }
 200   return (q);
 201 }
 202
 203
 204
 205 /* NAME
 206  *   filter_spatial_response_analytic                   Calculate filter by analytic inverse fourier
 207  *                              transform of filters's frequency
 208  *                              response
 209  *
 210  * SYNOPSIS
 211  *   y = filter_spatial_response_analytic (filt_type, x, bw, param)
 212  *   double y                   Filter's response in spatial domain
 213  *   int filt_type              Type of filter (definitions in ct.h)
 214  *   double x                   Spatial position to evaluate filter
 215  *   double bw                  Bandwidth of window
 216  *   double param               General parameter for various filters
 217  */
 218
 219 double
 220 filter_spatial_response_analytic (int filt_type, double x, double bw, double param)
 221 {
 222   double q, temp;
 223   double u = TWOPI * x;
 224   double w = bw / 2;
 225   double b = PI / bw;
 226   double b2 = TWOPI / bw;
 227
 228   switch (filt_type) {
 229   case FILTER_BANDLIMIT:
 230     q = bw * sinc(u * w, 1.0);
 231     break;
 232   case FILTER_TRIANGLE:
 233     temp = sinc (u * w, 1.0);
 234     q = bw * temp * temp;
 235     break;
 236   case FILTER_COSINE:
 237     q = sinc(b-u,w) + sinc(b+u,w);
 238     break;
 239   case FILTER_G_HAMMING:
 240     q = 2 * param * sin(u*w)/u + (1-param) *
 241       (sinc(b2-u, w) + sinc(b2+u, w));
 242     break;
 243   case FILTER_ABS_BANDLIMIT:
 244     q = 2 * integral_abscos (u, w);
 245     break;
 246   case FILTER_ABS_COSINE:
 247     q = integral_abscos(b-u,w) + integral_abscos(b+u,w);
 248     break;
 249   case FILTER_ABS_G_HAMMING:
 250     q = 2 * param * integral_abscos(u,w) +
 251       (1-param)*(integral_abscos(u-b2,w)+integral_abscos(u+b2,w));
 252     break;
 253   case FILTER_SHEPP:
 254     if (fabs (u) < 1E-7)
 255       q = 4. / (PI * bw * bw);
 256     else
 257       q = fabs ((2 / bw) * sin (u * w)) * sinc (u * w, 1.) * sinc (u * w, 1.);
 258     break;
 259   case FILTER_SINC:
 260     if (fabs (x) < bw / 2)
 261       q = 1.;
 262     else
 263       q = 0.;
 264     break;
 265   case FILTER_ABS_SINC:
 266   default:
 267     sys_error (ERR_WARNING,
 268                "Analytic filter type %d not implemented [filter_spatial_response_analytic]",
 269                filt_type);
 270     q = 0;
 271     break;
 272   }
 273
 274   return (q);
 275 }
 276
 277
 278 /* NAME
 279  *   sinc                       Return sin(x)/x function
 280  *
 281  * SYNOPSIS
 282  *   v = sinc (x, mult)
 283  *   double v                   sinc value
 284  *   double x, mult
 285  *
 286  * DESCRIPTION
 287  *   v = sin(x * mult) / x;
 288  */
 289
 290 double
 291 sinc (double x, double mult)
 292 {
 293   return (fabs(x) > F_EPSILON ? (sin (x * mult) / x) : 1.0);
 294 }
 295
 296
 297 /* NAME
 298  *   integral_abscos                    Returns integral of u*cos(u)
 299  *
 300  * SYNOPSIS
 301  *   q = integral_abscos (u, w)
 302  *   double q                   Integral value
 303  *   double u                   Integration variable
 304  *   double w                   Upper integration boundary
 305  *
 306  * DESCRIPTION
 307  *   Returns the value of integral of u*cos(u)*dV for V = 0 to w
 308  */
 309
 310 double
 311 integral_abscos (double u, double w)
 312 {
 313   if (fabs (u) > F_EPSILON)
 314     return (cos(u * w) - 1) / (u * u) + w / u * sin (u * w);
 315   else
 316     return (w * w / 2);
 317 }
 318
 319