3 This file is part of a program that implements a Software-Defined Radio.
5 Copyright (C) 2004 by Frank Brickle, AB2KT and Bob McGwier, N4HY
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.
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.
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
21 The authors can be reached by email at
29 The DTTS Microwave Society
43 int adaptive_filter_size,
45 LMSR lms = (LMSR) safealloc(1, sizeof(_lmsstate), "new_lmsr state");
48 lms->signal_size = CXBsize(lms->signal);
51 lms->mask = lms->size - 1;
52 lms->delay_line = newvec_REAL(lms->size, "lmsr delay");
53 lms->adaptation_rate = adaptation_rate;
54 lms->leakage = leakage;
55 lms->adaptive_filter_size = adaptive_filter_size;
56 lms->adaptive_filter = newvec_REAL(128, "lmsr filter");
57 lms->filter_type = filter_type;
58 lms->delay_line_ptr = 0;
66 delvec_REAL(lms->delay_line);
67 delvec_REAL(lms->adaptive_filter);
68 safefree((char *) lms);
72 // just to make the algorithm itself a little clearer,
73 // get the admin stuff out of the way
75 #define ssiz (lms->signal_size)
76 #define asiz (lms->adaptive_filter_size)
77 #define dptr (lms->delay_line_ptr)
78 #define rate (lms->adaptation_rate)
79 #define leak (lms->leakage)
81 #define ssig(n) (CXBreal(lms->signal,(n)))
83 #define dlay(n) (lms->delay_line[(n)])
85 #define afil(n) (lms->adaptive_filter[(n)])
86 #define wrap(n) (((n) + (lms->delay) + (lms->delay_line_ptr)) & (lms->mask))
87 #define bump(n) (((n) + (lms->mask)) & (lms->mask))
90 lmsr_adapt_i(LMSR lms) {
92 REAL sum_sq, scl1, scl2;
95 scl1 = 1.0 - rate * leak;
97 for (i = 0; i < ssiz; i++) {
103 for (j = 0; j < asiz; j++) {
105 sum_sq += sqr(dlay(k));
106 accum += afil(j) * dlay(k);
109 error = ssig(i) - accum;
112 scl2 = rate / (sum_sq + 1e-10);
114 for (j = 0; j < asiz; j++) {
116 afil(j) = afil(j) * scl1 + error * dlay(k);
124 lmsr_adapt_n(LMSR lms) {
126 REAL sum_sq, scl1, scl2;
129 scl1 = 1.0 - rate * leak;
131 for (i = 0; i < ssiz; i++) {
133 dlay(dptr) = ssig(i);
137 for (j = 0; j < asiz; j++) {
139 sum_sq += sqr(dlay(k));
140 accum += afil(j) * dlay(k);
143 error = ssig(i) - accum;
146 scl2 = rate / (sum_sq + 1e-10);
148 for (j = 0; j < asiz; j++) {
150 afil(j) = afil(j) * scl1 + error * dlay(k);
162 REAL adaptation_rate,
164 int adaptive_filter_size,
166 LMSR lms = (LMSR) safealloc(1, sizeof(_lmsstate), "new_lmsr state");
168 lms->signal = signal;
169 CXBsize(lms->signal);
172 lms->mask = lms->size - 1;
173 lms->delay_line = newvec_COMPLEX(lms->size, "lmsr delay");
174 lms->adaptation_rate = adaptation_rate;
175 lms->leakage = leakage;
176 lms->adaptive_filter_size = adaptive_filter_size;
177 lms->adaptive_filter = newvec_COMPLEX(128, "lmsr filter");
178 lms->filter_type = filter_type;
179 lms->delay_line_ptr = 0;
187 delvec_COMPLEX(lms->delay_line);
188 delvec_COMPLEX(lms->adaptive_filter);
189 safefree((char *) lms);
193 // just to make the algorithm itself a little clearer,
194 // get the admin stuff out of the way
196 #define ssiz (lms->signal_size)
197 #define asiz (lms->adaptive_filter_size)
198 #define dptr (lms->delay_line_ptr)
199 #define rate (lms->adaptation_rate)
200 #define leak (lms->leakage)
202 #define ssig(n) (CXBdata(lms->signal,(n)))
204 #define dlay(n) (lms->delay_line[(n)])
206 #define afil(n) (lms->adaptive_filter[(n)])
207 #define wrap(n) (((n) + (lms->delay) + (lms->delay_line_ptr)) & (lms->mask))
208 #define bump(n) (((n) + (lms->mask)) & (lms->mask))
211 lmsr_adapt_i(LMSR lms) {
213 REAL sum_sq, scl1, scl2;
214 COMPLEX accum, error;
215 scl1 = 1.0 - rate * leak;
216 for (i = 0; i < ssiz; i++) {
218 dlay(dptr) = ssig(i);
222 for (j = 0; j < asiz; j++) {
224 sum_sq += Csqrmag(dlay(k));
225 accum = Cadd(accum, Cmul(Conjg(afil(j)), dlay(k)));
228 error = Csub(ssig(i), accum);
231 scl2 = rate / (sum_sq + 1e-10);
232 error = Cscl(Conjg(error),scl2);
233 for (j = 0; j < asiz; j++) {
235 afil(j) = Cadd(Cscl(afil(j), scl1),Cmul(error, dlay(k)));
243 lmsr_adapt_n(LMSR lms) {
245 REAL sum_sq, scl1, scl2;
246 COMPLEX accum, error;
247 scl1 = 1.0 - rate * leak;
248 for (i = 0; i < ssiz; i++) {
250 dlay(dptr) = ssig(i);
254 for (j = 0; j < asiz; j++) {
256 sum_sq += Csqrmag(dlay(k));
257 accum = Cadd(accum, Cmul(Conjg(afil(j)), dlay(k)));
260 error = Csub(ssig(i), accum);
263 scl2 = rate / (sum_sq + 1e-10);
264 error = Cscl(Conjg(error),scl2);
265 for (j = 0; j < asiz; j++) {
267 afil(j) = Cadd(Cscl(afil(j), scl1),Cmul(error, dlay(k)));
277 lmsr_adapt(LMSR lms) {
278 switch(lms->filter_type) {
282 case LMADF_INTERFERENCE: