]> git.rkrishnan.org Git - dttsp.git/blob - jDttSP/ovsv.c
Upgraded windows alternatives
[dttsp.git] / 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