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 void
  96 reset_OvSv(FiltOvSv pflt) {
  97   memset((char *) pflt->zrvec, 0, pflt->fftlen * sizeof(COMPLEX));
  98 }
  99
 100 /*------------------------------------------------------------*/
 101 /* info: */
 102 /* NB strategy. This is the address we pass to newCXB as
 103    the place to read samples into. It's the right half of
 104    the true buffer. Old samples get slid from here to
 105    left half after each go-around. */
 106 COMPLEX *
 107 FiltOvSv_initpoint(FiltOvSv pflt) { return &(pflt->zrvec[pflt->buflen]); }
 108
 109 /* how many to put there */
 110 int
 111 FiltOvSv_initsize(FiltOvSv pflt) { return (pflt->fftlen - pflt->buflen); }
 112
 113 /* where to put next batch of samples to filter */
 114 COMPLEX *
 115 FiltOvSv_fetchpoint(FiltOvSv pflt) { return &(pflt->zrvec[pflt->buflen]); }
 116
 117 /* how many samples to put there */
 118
 119 int
 120 FiltOvSv_fetchsize(FiltOvSv pflt) { return (pflt->fftlen - pflt->buflen); }
 121
 122 /* where samples should be taken from after filtering */
 123 #ifdef LHS
 124 COMPLEX *
 125 FiltOvSv_storepoint(FiltOvSv plft) { return ((pflt->zovec) + buflen); }
 126 #else
 127 COMPLEX *
 128 FiltOvSv_storepoint(FiltOvSv pflt) { return ((pflt->zovec)); }
 129 #endif
 130
 131 /* alternating parity fetch/store */
 132
 133 COMPLEX *
 134 FiltOvSv_fetchpt_par(FiltOvSv pflt, int parity) {
 135   if (parity & 01)
 136     return pflt->zovec + pflt->buflen;
 137   else
 138     return pflt->zovec;
 139 }
 140
 141 COMPLEX *
 142 FiltOvSv_storept_par(FiltOvSv pflt, int parity) {
 143   if (parity & 01)
 144     return pflt->zovec;
 145   else
 146     return pflt->zovec + pflt->buflen;
 147 }
 148
 149 /* how many samples to take */
 150 /* NB strategy. This is the number of good samples in the
 151    left half of the true buffer. Samples in right half
 152    are circular artifacts and are ignored. */
 153 int
 154 FiltOvSv_storesize(FiltOvSv pflt) { return (pflt->fftlen - pflt->buflen); }
 155
 156 /*------------------------------------------------------------*/
 157 /* create a new overlap/save filter from complex coefficients */
 158
 159 FiltOvSv
 160 newFiltOvSv(COMPLEX *coefs, int ncoef, int pbits) {
 161   int buflen, fftlen;
 162   FiltOvSv p;
 163   fftw_plan pfwd,pinv;
 164   COMPLEX *zrvec, *zfvec, *zivec, *zovec;
 165   p = (FiltOvSv) safealloc(1, sizeof(filt_ov_sv), "new overlap/save filter");
 166   buflen = nblock2(ncoef-1), fftlen = 2 * buflen;
 167   zrvec = newvec_COMPLEX(fftlen, "raw signal vec in newFiltOvSv");
 168   zfvec = newvec_COMPLEX(fftlen, "filter z vec in newFiltOvSv");
 169   zivec = newvec_COMPLEX(fftlen, "signal in z vec in newFiltOvSv");
 170   zovec = newvec_COMPLEX(fftlen, "signal out z vec in newFiltOvSv");
 171
 172   /* prepare frequency response from filter coefs */
 173   {
 174     int i;
 175     COMPLEX *zcvec;
 176     fftw_plan ptmp;
 177
 178     zcvec = newvec_COMPLEX(fftlen, "temp filter z vec in newFiltOvSv");
 179     ptmp = fftw_create_plan(fftlen, FFTW_FORWARD, pbits);
 180
 181 #ifdef LHS
 182     for (i = 0; i < ncoef; i++) zcvec[i] = coefs[i];
 183 #else
 184     for (i = 0; i < ncoef; i++) zcvec[fftlen - ncoef - 1 + i] = coefs[i];
 185 #endif
 186
 187     fftw_one(ptmp, (fftw_complex *) zcvec, (fftw_complex *) zfvec);
 188     fftw_destroy_plan(ptmp);
 189     delvec_COMPLEX(zcvec);
 190   }
 191
 192   /* prepare transforms for signal */
 193   pfwd = fftw_create_plan(fftlen, FFTW_FORWARD, pbits);
 194   pinv = fftw_create_plan(fftlen, FFTW_BACKWARD, pbits);
 195   /* stuff values */
 196   p->buflen = buflen;
 197   p->fftlen = fftlen;
 198   p->zfvec = zfvec;
 199   p->zivec = zivec;
 200   p->zovec = zovec;
 201   p->zrvec = zrvec;
 202   p->pfwd = pfwd;
 203   p->pinv = pinv;
 204   p->scale = 1.0 / fftlen;
 205
 206   return p;
 207 }
 208
 209 /* deep-six the filter */
 210 void
 211 delFiltOvSv(FiltOvSv p) {
 212   if (p) {
 213     delvec_COMPLEX(p->zfvec);
 214     delvec_COMPLEX(p->zivec);
 215     delvec_COMPLEX(p->zovec);
 216     delvec_COMPLEX(p->zrvec);
 217     fftw_destroy_plan(p->pfwd);
 218     fftw_destroy_plan(p->pinv);
 219     safefree((char *) p);
 220   }
 221 }
 222
 223 /*------------------------------------------------------------*/
 224
 225