Qucs-core  0.0.19
receiver.cpp
Go to the documentation of this file.
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