Qucs-core
0.0.19
|
00001 /* 00002 * receiver.cpp - receiver transformation class implementation 00003 * 00004 * Copyright (C) 2009 Dirk Schaefer <schad@5pm.de> 00005 * Copyright (C) 2009 Stefan Jahn <stefan@lkcc.org> 00006 * 00007 * This is free software; you can redistribute it and/or modify 00008 * it under the terms of the GNU General Public License as published by 00009 * the Free Software Foundation; either version 2, or (at your option) 00010 * any later version. 00011 * 00012 * This software is distributed in the hope that it will be useful, 00013 * but WITHOUT ANY WARRANTY; without even the implied warranty of 00014 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 00015 * GNU General Public License for more details. 00016 * 00017 * You should have received a copy of the GNU General Public License 00018 * along with this package; see the file COPYING. If not, write to 00019 * the Free Software Foundation, Inc., 51 Franklin Street - Fifth Floor, 00020 * Boston, MA 02110-1301, USA. 00021 * 00022 * $Id$ 00023 * 00024 */ 00025 00026 #if HAVE_CONFIG_H 00027 # include <config.h> 00028 #endif 00029 00030 #include <stdio.h> 00031 #include <stdlib.h> 00032 #include <string.h> 00033 #include <cmath> 00034 #include <cstdint> 00035 00036 #include "consts.h" 00037 #include "object.h" 00038 #include "complex.h" 00039 #include "vector.h" 00040 #include "spline.h" 00041 #include "interpolator.h" 00042 #include "fourier.h" 00043 #include "receiver.h" 00044 00045 namespace qucs { 00046 00047 /* The function returns a power-of-two value which is equal or larger 00048 than the given value. The maximum returned value is 2^30. */ 00049 int32_t emi::nearestbin32 (int x) { 00050 int32_t boundary = 1 << 30; 00051 int32_t current = 1; 00052 if (x >= boundary) return boundary; 00053 while (current < x) current <<= 1; 00054 return current; 00055 } 00056 00057 /* Ideal filter construction for given center frequency, bandwidth and 00058 the frequency for which the filter is evaluated. */ 00059 nr_double_t emi::f_ideal (nr_double_t fc, nr_double_t bw, nr_double_t f) { 00060 nr_double_t lo = fc - bw / 2; 00061 nr_double_t hi = fc + bw / 2; 00062 if (f >= lo && f < hi) 00063 return 1.0; 00064 return 0.0; 00065 } 00066 00067 /* Construction of a bandpass filter of 2nd order for given center 00068 frequency, bandwidth and the frequency for which the filter is 00069 evaluated. */ 00070 nr_double_t emi::f_2ndorder (nr_double_t fc, nr_double_t bw, nr_double_t f) { 00071 nr_double_t q = fc / bw; 00072 nr_complex_t p = nr_complex_t (0, f / fc); 00073 nr_complex_t w = p / q / (1.0 + p / q + p * p); 00074 return norm (w); 00075 } 00076 00077 /* Construction of a gaussian filter for given center frequency, 00078 bandwidth and the frequency for which the filter is evaluated. */ 00079 nr_double_t emi::f_gauss (nr_double_t fc, nr_double_t bw, nr_double_t f) { 00080 nr_double_t a = log (0.5) / bw / bw; 00081 nr_double_t s = f - fc; 00082 return exp (a * s * s); 00083 } 00084 00085 /* The function computes a EMI receiver spectrum based on the given 00086 waveform in the time domain. The number of points in the waveform 00087 is required to be a power of two. Also the samples are supposed 00088 to be equidistant. */ 00089 vector * emi::receiver (nr_double_t * ida, nr_double_t duration, int ilength) { 00090 00091 int i, n, points; 00092 nr_double_t fres; 00093 vector * ed = new vector (); 00094 00095 points = ilength; 00096 00097 /* ilength must be a power of 2 - write wrapper later on */ 00098 fourier::_fft_1d (ida, ilength, 1); /* 1 = forward fft */ 00099 00100 /* start at first AC point (0 as DC point remains untouched) 00101 additionally only half of the FFT result required */ 00102 for (i = 2; i < points; i++) { 00103 ida[i] /= points / 2; 00104 } 00105 00106 /* calculate frequency step */ 00107 fres = 1.0 / duration; 00108 00109 /* generate data vector; inplace calculation of magnitudes */ 00110 nr_double_t * d = ida; 00111 for (n = 0, i = 0; i < points / 2; i++, n += 2){ 00112 /* abs value of complex number */ 00113 d[i] = xhypot (ida[n], ida[n + 1]); 00114 /* vector contains complex values; thus every second value */ 00115 } 00116 00117 points /= 2; 00118 00119 /* define EMI settings */ 00120 struct settings settings[] = { 00121 { 200, 150e3, 200, 200 }, 00122 { 150e3, 30e6, 9e3, 9e3 }, 00123 { 30e6, 1e9, 120e3, 120e3 }, 00124 { 0, 0, 0, 0} 00125 }; 00126 00127 /* define EMI noise floor */ 00128 nr_double_t noise = std::pow (10.0, (-100.0 / 40.0)) * 1e-6; 00129 00130 /* generate resulting data & frequency vector */ 00131 nr_double_t fcur, dcur; 00132 int ei = 0; 00133 00134 /* go through each EMI setting */ 00135 for (i = 0; settings[i].bandwidth != 0; i++ ) { 00136 00137 nr_double_t bw = settings[i].bandwidth; 00138 nr_double_t fstart = settings[i].start; 00139 nr_double_t fstop = settings[i].stop; 00140 nr_double_t fstep = settings[i].stepsize; 00141 00142 /* go through frequencies */ 00143 for (fcur = fstart; fcur <= fstop; fcur += fstep) { 00144 00145 /* calculate upper and lower frequency bounds */ 00146 nr_double_t lo = fcur - bw / 2; 00147 nr_double_t hi = fcur + bw / 2; 00148 if (hi < fres) continue; 00149 00150 /* calculate indices covering current bandwidth */ 00151 int il = std::floor (lo / fres); 00152 int ir = std::floor (hi / fres); 00153 00154 /* right index (ri) greater 0 and left index less than points -> 00155 at least part of data is within bandwidth indices */ 00156 if (ir >= 0 && il < points - 1) { 00157 /* adjust indices to reside in the data array */ 00158 if (il < 0) il = 0; 00159 if (ir > points - 1) ir = points - 1; 00160 00161 /* sum-up the values within the bandwidth */ 00162 dcur = 0; 00163 for (int j = 0; j < ir - il; j++){ 00164 nr_double_t f = fres * (il + j); 00165 dcur += f_2ndorder (fcur, bw, f) * d[il + j]; 00166 } 00167 00168 /* add noise to the result */ 00169 dcur += noise * sqrt (bw); 00170 00171 /* save result */ 00172 ed->add (nr_complex_t (dcur, fcur)); 00173 ei++; 00174 } 00175 } 00176 } 00177 00178 /* returning values of function */ 00179 return ed; 00180 } 00181 00182 /* This is a wrapper for the basic EMI rceiver functionality. It 00183 takes an arbitraty waveform in the time domain and interpolates it 00184 such, that its length results in a power of two elements. */ 00185 vector * emi::receiver (vector * da, vector * dt, int len) { 00186 00187 int i, nlen, olen = da->getSize (); 00188 00189 // don't allow less points than the actual length 00190 if (len < da->getSize ()) len = da->getSize (); 00191 00192 // find a power-of-two length 00193 nlen = emi::nearestbin32 (len); 00194 00195 nr_double_t tstart = real (dt->get (0)); 00196 nr_double_t tstop = real (dt->get (olen - 1)); 00197 nr_double_t duration = tstop - tstart; 00198 00199 /* please note: interpolation is always performed in order to ensure 00200 equidistant samples */ 00201 00202 // create interpolator (use cubic splines) 00203 interpolator * inter = new interpolator (); 00204 inter->rvectors (da, dt); 00205 inter->prepare (INTERPOL_CUBIC, REPEAT_NO, DATA_RECTANGULAR); 00206 00207 // adjust the time domain vector using interpolation 00208 nr_double_t * ida = new nr_double_t[2 * nlen]; 00209 nr_double_t tstep = duration / (nlen - 1); 00210 for (i = 0; i < nlen; i++) { 00211 nr_double_t t = i * tstep + tstart; 00212 ida[2 * i + 0] = inter->rinterpolate (t); 00213 ida[2 * i + 1] = 0; 00214 } 00215 00216 // destroy interpolator 00217 delete inter; 00218 00219 // run actual EMI receiver computations 00220 vector * res = receiver (ida, duration, nlen); 00221 00222 // delete intermediate data 00223 delete[] ida; 00224 00225 return res; 00226 } 00227 00228 } // namespace qucs