Branch data Line data Source code
1 : : /*
2 : : * receiver.cpp - receiver transformation class implementation
3 : : *
4 : : * Copyright (C) 2009 Dirk Schaefer <schad@5pm.de>
5 : : * Copyright (C) 2009 Stefan Jahn <stefan@lkcc.org>
6 : : *
7 : : * This 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, or (at your option)
10 : : * any later version.
11 : : *
12 : : * This software 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 package; see the file COPYING. If not, write to
19 : : * the Free Software Foundation, Inc., 51 Franklin Street - Fifth Floor,
20 : : * Boston, MA 02110-1301, USA.
21 : : *
22 : : * $Id$
23 : : *
24 : : */
25 : :
26 : : #if HAVE_CONFIG_H
27 : : # include <config.h>
28 : : #endif
29 : :
30 : : #include <stdio.h>
31 : : #include <stdlib.h>
32 : : #include <string.h>
33 : : #include <cmath>
34 : : #include <cstdint>
35 : :
36 : : #include "consts.h"
37 : : #include "object.h"
38 : : #include "complex.h"
39 : : #include "vector.h"
40 : : #include "spline.h"
41 : : #include "interpolator.h"
42 : : #include "fourier.h"
43 : : #include "receiver.h"
44 : :
45 : : namespace qucs {
46 : :
47 : : /* The function returns a power-of-two value which is equal or larger
48 : : than the given value. The maximum returned value is 2^30. */
49 : 0 : int32_t emi::nearestbin32 (int x) {
50 : 0 : int32_t boundary = 1 << 30;
51 : 0 : int32_t current = 1;
52 [ # # ]: 0 : if (x >= boundary) return boundary;
53 [ # # ]: 0 : while (current < x) current <<= 1;
54 : 0 : return current;
55 : : }
56 : :
57 : : /* Ideal filter construction for given center frequency, bandwidth and
58 : : the frequency for which the filter is evaluated. */
59 : 0 : nr_double_t emi::f_ideal (nr_double_t fc, nr_double_t bw, nr_double_t f) {
60 : 0 : nr_double_t lo = fc - bw / 2;
61 : 0 : nr_double_t hi = fc + bw / 2;
62 [ # # ][ # # ]: 0 : if (f >= lo && f < hi)
63 : 0 : return 1.0;
64 : 0 : return 0.0;
65 : : }
66 : :
67 : : /* Construction of a bandpass filter of 2nd order for given center
68 : : frequency, bandwidth and the frequency for which the filter is
69 : : evaluated. */
70 : 0 : nr_double_t emi::f_2ndorder (nr_double_t fc, nr_double_t bw, nr_double_t f) {
71 : 0 : nr_double_t q = fc / bw;
72 : 0 : nr_complex_t p = nr_complex_t (0, f / fc);
73 [ # # ][ # # ]: 0 : nr_complex_t w = p / q / (1.0 + p / q + p * p);
[ # # ]
74 [ # # ]: 0 : return norm (w);
75 : : }
76 : :
77 : : /* Construction of a gaussian filter for given center frequency,
78 : : bandwidth and the frequency for which the filter is evaluated. */
79 : 0 : nr_double_t emi::f_gauss (nr_double_t fc, nr_double_t bw, nr_double_t f) {
80 : 0 : nr_double_t a = log (0.5) / bw / bw;
81 : 0 : nr_double_t s = f - fc;
82 : 0 : return exp (a * s * s);
83 : : }
84 : :
85 : : /* The function computes a EMI receiver spectrum based on the given
86 : : waveform in the time domain. The number of points in the waveform
87 : : is required to be a power of two. Also the samples are supposed
88 : : to be equidistant. */
89 : 0 : vector * emi::receiver (nr_double_t * ida, nr_double_t duration, int ilength) {
90 : :
91 : : int i, n, points;
92 : : nr_double_t fres;
93 [ # # ][ # # ]: 0 : vector * ed = new vector ();
94 : :
95 : 0 : points = ilength;
96 : :
97 : : /* ilength must be a power of 2 - write wrapper later on */
98 [ # # ]: 0 : fourier::_fft_1d (ida, ilength, 1); /* 1 = forward fft */
99 : :
100 : : /* start at first AC point (0 as DC point remains untouched)
101 : : additionally only half of the FFT result required */
102 [ # # ]: 0 : for (i = 2; i < points; i++) {
103 : 0 : ida[i] /= points / 2;
104 : : }
105 : :
106 : : /* calculate frequency step */
107 : 0 : fres = 1.0 / duration;
108 : :
109 : : /* generate data vector; inplace calculation of magnitudes */
110 : 0 : nr_double_t * d = ida;
111 [ # # ]: 0 : for (n = 0, i = 0; i < points / 2; i++, n += 2){
112 : : /* abs value of complex number */
113 [ # # ]: 0 : d[i] = xhypot (ida[n], ida[n + 1]);
114 : : /* vector contains complex values; thus every second value */
115 : : }
116 : :
117 : 0 : points /= 2;
118 : :
119 : : /* define EMI settings */
120 : : struct settings settings[] = {
121 : : { 200, 150e3, 200, 200 },
122 : : { 150e3, 30e6, 9e3, 9e3 },
123 : : { 30e6, 1e9, 120e3, 120e3 },
124 : : { 0, 0, 0, 0}
125 : 0 : };
126 : :
127 : : /* define EMI noise floor */
128 : 0 : nr_double_t noise = std::pow (10.0, (-100.0 / 40.0)) * 1e-6;
129 : :
130 : : /* generate resulting data & frequency vector */
131 : : nr_double_t fcur, dcur;
132 : 0 : int ei = 0;
133 : :
134 : : /* go through each EMI setting */
135 [ # # ]: 0 : for (i = 0; settings[i].bandwidth != 0; i++ ) {
136 : :
137 : 0 : nr_double_t bw = settings[i].bandwidth;
138 : 0 : nr_double_t fstart = settings[i].start;
139 : 0 : nr_double_t fstop = settings[i].stop;
140 : 0 : nr_double_t fstep = settings[i].stepsize;
141 : :
142 : : /* go through frequencies */
143 [ # # ]: 0 : for (fcur = fstart; fcur <= fstop; fcur += fstep) {
144 : :
145 : : /* calculate upper and lower frequency bounds */
146 : 0 : nr_double_t lo = fcur - bw / 2;
147 : 0 : nr_double_t hi = fcur + bw / 2;
148 [ # # ]: 0 : if (hi < fres) continue;
149 : :
150 : : /* calculate indices covering current bandwidth */
151 : 0 : int il = std::floor (lo / fres);
152 : 0 : int ir = std::floor (hi / fres);
153 : :
154 : : /* right index (ri) greater 0 and left index less than points ->
155 : : at least part of data is within bandwidth indices */
156 [ # # ][ # # ]: 0 : if (ir >= 0 && il < points - 1) {
157 : : /* adjust indices to reside in the data array */
158 [ # # ]: 0 : if (il < 0) il = 0;
159 [ # # ]: 0 : if (ir > points - 1) ir = points - 1;
160 : :
161 : : /* sum-up the values within the bandwidth */
162 : 0 : dcur = 0;
163 [ # # ]: 0 : for (int j = 0; j < ir - il; j++){
164 : 0 : nr_double_t f = fres * (il + j);
165 [ # # ]: 0 : dcur += f_2ndorder (fcur, bw, f) * d[il + j];
166 : : }
167 : :
168 : : /* add noise to the result */
169 [ # # ]: 0 : dcur += noise * sqrt (bw);
170 : :
171 : : /* save result */
172 [ # # ]: 0 : ed->add (nr_complex_t (dcur, fcur));
173 : 0 : ei++;
174 : : }
175 : : }
176 : : }
177 : :
178 : : /* returning values of function */
179 : 0 : return ed;
180 : : }
181 : :
182 : : /* This is a wrapper for the basic EMI rceiver functionality. It
183 : : takes an arbitraty waveform in the time domain and interpolates it
184 : : such, that its length results in a power of two elements. */
185 : 0 : vector * emi::receiver (vector * da, vector * dt, int len) {
186 : :
187 : 0 : int i, nlen, olen = da->getSize ();
188 : :
189 : : // don't allow less points than the actual length
190 [ # # ]: 0 : if (len < da->getSize ()) len = da->getSize ();
191 : :
192 : : // find a power-of-two length
193 : 0 : nlen = emi::nearestbin32 (len);
194 : :
195 : 0 : nr_double_t tstart = real (dt->get (0));
196 : 0 : nr_double_t tstop = real (dt->get (olen - 1));
197 : 0 : nr_double_t duration = tstop - tstart;
198 : :
199 : : /* please note: interpolation is always performed in order to ensure
200 : : equidistant samples */
201 : :
202 : : // create interpolator (use cubic splines)
203 [ # # ]: 0 : interpolator * inter = new interpolator ();
204 : 0 : inter->rvectors (da, dt);
205 : 0 : inter->prepare (INTERPOL_CUBIC, REPEAT_NO, DATA_RECTANGULAR);
206 : :
207 : : // adjust the time domain vector using interpolation
208 : 0 : nr_double_t * ida = new nr_double_t[2 * nlen];
209 : 0 : nr_double_t tstep = duration / (nlen - 1);
210 [ # # ]: 0 : for (i = 0; i < nlen; i++) {
211 : 0 : nr_double_t t = i * tstep + tstart;
212 : 0 : ida[2 * i + 0] = inter->rinterpolate (t);
213 : 0 : ida[2 * i + 1] = 0;
214 : : }
215 : :
216 : : // destroy interpolator
217 [ # # ]: 0 : delete inter;
218 : :
219 : : // run actual EMI receiver computations
220 : 0 : vector * res = receiver (ida, duration, nlen);
221 : :
222 : : // delete intermediate data
223 [ # # ]: 0 : delete[] ida;
224 : :
225 : 0 : return res;
226 : : }
227 : :
228 : : } // namespace qucs
|
