jDttSP/ovsv.c

   1 /* ovsv.c
   2
   3 This file is part of a program that implements a Software-Defined Radio.
   4
   5 Copyright (C) 2004 by Frank Brickle, AB2KT and Bob McGwier, N4HY
   6
   7 This program is free software; you can redistribute it and/or modify
   8 it under the terms of the GNU General Public License as published by
   9 the Free Software Foundation; either version 2 of the License, or
  10 (at your option) any later version.
  11
  12 This program is distributed in the hope that it will be useful,
  13 but WITHOUT ANY WARRANTY; without even the implied warranty of
  14 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  15 GNU General Public License for more details.
  16
  17 You should have received a copy of the GNU General Public License
  18 along with this program; if not, write to the Free Software
  19 Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
  20
  21 The authors can be reached by email at
  22
  23 ab2kt@arrl.net
  24 or
  25 rwmcgwier@comcast.net
  26
  27 or by paper mail at
  28
  29 The DTTS Microwave Society
  30 6 Kathleen Place
  31 Bridgewater, NJ 08807
  32 */
  33
  34 #include <common.h>
  35
  36 /*------------------------------------------------------------*/
  37 /* run the filter */
  38
  39 void
  40 filter_OvSv(FiltOvSv pflt) {
  41   int i,
  42       j,
  43       m = pflt->fftlen,
  44       n = pflt->buflen;
  45   COMPLEX *zfvec = pflt->zfvec,
  46           *zivec = pflt->zivec,
  47           *zovec = pflt->zovec,
  48           *zrvec = pflt->zrvec;
  49   REAL scl = pflt->scale;
  50
  51   /* input sig -> z */
  52   fftw_one(pflt->pfwd, (fftw_complex *) zrvec, (fftw_complex *) zivec);
  53
  54   /* convolve in z */
  55   for (i = 0; i < m; i++) zivec[i] = Cmul(zivec[i], zfvec[i]);
  56
  57   /* z convolved sig -> time output sig */
  58   fftw_one(pflt->pinv, (fftw_complex *) zivec, (fftw_complex *) zovec);
  59
  60   /* scale */
  61   for (i = 0; i < n; i++) zovec[i].re *= scl, zovec[i].im *= scl;
  62
  63   /* prepare input sig vec for next fill */
  64   for (i = 0, j = n; i < n; i++, j++) zrvec[i] = zrvec[j];
  65 }
  66
  67 void
  68 filter_OvSv_par(FiltOvSv pflt) {
  69   int i,
  70       j,
  71       m = pflt->fftlen,
  72       n = pflt->buflen;
  73   COMPLEX *zfvec = pflt->zfvec,
  74           *zivec = pflt->zivec,
  75           *zovec = pflt->zovec,
  76           *zrvec = pflt->zrvec;
  77   REAL scl = pflt->scale;
  78
  79   /* input sig -> z */
  80   fftw_one(pflt->pfwd,(fftw_complex *) zrvec, (fftw_complex *) zivec);
  81
  82   /* convolve in z */
  83   for (i = 0; i < m; i++) zivec[i] = Cmul(zivec[i], zfvec[i]);
  84
  85   /* z convolved sig -> time output sig */
  86   fftw_one(pflt->pinv, (fftw_complex *) zivec, (fftw_complex *) zovec);
  87
  88   /* scale */
  89   for (i = 0; i < n; i++) zovec[i].re *= scl, zovec[i].im *= scl;
  90
  91   /* prepare input sig vec for next fill */
  92   for (i = 0, j = n; i < n; i++, j++) zrvec[i] = zrvec[j];
  93 }
  94
  95 /*------------------------------------------------------------*/
  96 /* info: */
  97 /* NB strategy. This is the address we pass to newCXB as
  98    the place to read samples into. It's the right half of
  99    the true buffer. Old samples get slid from here to
 100    left half after each go-around. */
 101 COMPLEX *
 102 FiltOvSv_initpoint(FiltOvSv pflt) { return &(pflt->zrvec[pflt->buflen]); }
 103
 104 /* how many to put there */
 105 int
 106 FiltOvSv_initsize(FiltOvSv pflt) { return (pflt->fftlen - pflt->buflen); }
 107
 108 /* where to put next batch of samples to filter */
 109 COMPLEX *
 110 FiltOvSv_fetchpoint(FiltOvSv pflt) { return &(pflt->zrvec[pflt->buflen]); }
 111
 112 /* how many samples to put there */
 113
 114 int
 115 FiltOvSv_fetchsize(FiltOvSv pflt) { return (pflt->fftlen - pflt->buflen); }
 116
 117 /* where samples should be taken from after filtering */
 118 #ifdef LHS
 119 COMPLEX *
 120 FiltOvSv_storepoint(FiltOvSv plft) { return ((pflt->zovec) + buflen); }
 121 #else
 122 COMPLEX *
 123 FiltOvSv_storepoint(FiltOvSv pflt) { return ((pflt->zovec)); }
 124 #endif
 125
 126 /* alternating parity fetch/store */
 127
 128 COMPLEX *
 129 FiltOvSv_fetchpt_par(FiltOvSv pflt, int parity) {
 130   if (parity & 01)
 131     return pflt->zovec + pflt->buflen;
 132   else
 133     return pflt->zovec;
 134 }
 135
 136 COMPLEX *
 137 FiltOvSv_storept_par(FiltOvSv pflt, int parity) {
 138   if (parity & 01)
 139     return pflt->zovec;
 140   else
 141     return pflt->zovec + pflt->buflen;
 142 }
 143
 144 /* how many samples to take */
 145 /* NB strategy. This is the number of good samples in the
 146    left half of the true buffer. Samples in right half
 147    are circular artifacts and are ignored. */
 148 int
 149 FiltOvSv_storesize(FiltOvSv pflt) { return (pflt->fftlen - pflt->buflen); }
 150
 151 /*------------------------------------------------------------*/
 152 /* create a new overlap/save filter from complex coefficients */
 153
 154 FiltOvSv
 155 newFiltOvSv(COMPLEX *coefs, int ncoef, int pbits) {
 156   int buflen, fftlen;
 157   FiltOvSv p;
 158   fftw_plan pfwd,pinv;
 159   COMPLEX *zrvec, *zfvec, *zivec, *zovec;
 160   p = (FiltOvSv) safealloc(1, sizeof(filt_ov_sv), "new overlap/save filter");
 161   buflen = nblock2(ncoef-1), fftlen = 2 * buflen;
 162   zrvec = newvec_COMPLEX(fftlen, "raw signal vec in newFiltOvSv");
 163   zfvec = newvec_COMPLEX(fftlen, "filter z vec in newFiltOvSv");
 164   zivec = newvec_COMPLEX(fftlen, "signal in z vec in newFiltOvSv");
 165   zovec = newvec_COMPLEX(fftlen, "signal out z vec in newFiltOvSv");
 166
 167   /* prepare frequency response from filter coefs */
 168   {
 169     int i;
 170     COMPLEX *zcvec;
 171     fftw_plan ptmp;
 172
 173     zcvec = newvec_COMPLEX(fftlen, "temp filter z vec in newFiltOvSv");
 174     ptmp = fftw_create_plan(fftlen, FFTW_FORWARD, pbits);
 175
 176 #ifdef LHS
 177     for (i = 0; i < ncoef; i++) zcvec[i] = coefs[i];
 178 #else
 179     for (i = 0; i < ncoef; i++) zcvec[fftlen - ncoef - 1 + i] = coefs[i];
 180 #endif
 181
 182     fftw_one(ptmp, (fftw_complex *) zcvec, (fftw_complex *) zfvec);
 183     fftw_destroy_plan(ptmp);
 184     delvec_COMPLEX(zcvec);
 185   }
 186
 187   /* prepare transforms for signal */
 188   pfwd = fftw_create_plan(fftlen, FFTW_FORWARD, pbits);
 189   pinv = fftw_create_plan(fftlen, FFTW_BACKWARD, pbits);
 190   /* stuff values */
 191   p->buflen = buflen;
 192   p->fftlen = fftlen;
 193   p->zfvec = zfvec;
 194   p->zivec = zivec;
 195   p->zovec = zovec;
 196   p->zrvec = zrvec;
 197   p->pfwd = pfwd;
 198   p->pinv = pinv;
 199   p->scale = 1.0 / fftlen;
 200
 201   return p;
 202 }
 203
 204 /* deep-six the filter */
 205 void
 206 delFiltOvSv(FiltOvSv p) {
 207   if (p) {
 208     delvec_COMPLEX(p->zfvec);
 209     delvec_COMPLEX(p->zivec);
 210     delvec_COMPLEX(p->zovec);
 211     delvec_COMPLEX(p->zrvec);
 212     fftw_destroy_plan(p->pfwd);
 213     fftw_destroy_plan(p->pinv);
 214     safefree((char *) p);
 215   }
 216 }
 217
 218 /*------------------------------------------------------------*/
 219
 220