Qucs-core  0.0.19
fourier.cpp
Go to the documentation of this file.
00001 /*
00002  * fourier.cpp - fourier transformation class implementation
00003  *
00004  * Copyright (C) 2005, 2006, 2009 Stefan Jahn <stefan@lkcc.org>
00005  *
00006  * This is free software; you can redistribute it and/or modify
00007  * it under the terms of the GNU General Public License as published by
00008  * the Free Software Foundation; either version 2, or (at your option)
00009  * any later version.
00010  *
00011  * This software is distributed in the hope that it will be useful,
00012  * but WITHOUT ANY WARRANTY; without even the implied warranty of
00013  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00014  * GNU General Public License for more details.
00015  *
00016  * You should have received a copy of the GNU General Public License
00017  * along with this package; see the file COPYING.  If not, write to
00018  * the Free Software Foundation, Inc., 51 Franklin Street - Fifth Floor,
00019  * Boston, MA 02110-1301, USA.
00020  *
00021  * $Id$
00022  *
00023  */
00024 
00025 #if HAVE_CONFIG_H
00026 # include <config.h>
00027 #endif
00028 
00029 #include <stdio.h>
00030 #include <stdlib.h>
00031 #include <string.h>
00032 #include <cmath>
00033 
00034 #include "consts.h"
00035 #include "object.h"
00036 #include "complex.h"
00037 #include "vector.h"
00038 #include "fourier.h"
00039 
00040 // Helper macro.
00041 #define Swap(a,b) { nr_double_t t; t = a; a = b; b = t; }
00042 
00043 namespace qucs {
00044 
00045 using namespace fourier;
00046 
00047 /* The function performs a 1-dimensional fast fourier transformation.
00048    Each data item is meant to be defined in equidistant steps.  The
00049    number of data items needs to be of binary size, e.g. 64, 128. */
00050 void fourier::_fft_1d (nr_double_t * data, int len, int isign) {
00051 
00052   // bit reversal method
00053   int i, j, m, n;
00054   n = 2 * len;
00055   j = 0;
00056   for (i = 0; i < n; i += 2) {
00057     if (j > i) {                   // was index already swapped ?
00058       Swap (data[j], data[i]);     // swap real part
00059       Swap (data[j+1], data[i+1]); // swap imaginary part
00060     }
00061     m = len;
00062     while (m >= 2 && j >= m) {     // calculate next swap index
00063       j -= m;
00064       m >>= 1;
00065     }
00066     j += m;
00067   }
00068 
00069   // Danielson-Lanzcos algorithm
00070   int mmax, istep;
00071   nr_double_t wt, th, wr, wi, wpr, wpi, ti, tr;
00072   mmax = 2;
00073   while (n > mmax) {
00074     istep = mmax << 1;
00075     th = isign * (2 * pi / mmax);
00076     wt = sin (0.5 * th);
00077     wpr = -2.0 * wt * wt;
00078     wpi = sin (th);
00079     wr = 1.0;
00080     wi = 0.0;
00081     for (m = 1; m < mmax; m += 2) {
00082       for (i = m; i <= n; i += istep) {
00083         j = i + mmax;
00084         tr = wr * data[j] - wi * data[j-1];
00085         ti = wr * data[j-1] + wi * data[j];
00086         data[j] = data[i] - tr;
00087         data[j-1] = data[i-1] - ti;
00088         data[i] += tr;
00089         data[i-1] += ti;
00090       }
00091       wr = (wt = wr) * wpr - wi * wpi + wr;
00092       wi = wi * wpr + wt * wpi + wi;
00093     }
00094     mmax = istep;
00095   }
00096 }
00097 
00098 /* The function transforms two real vectors using a single fast
00099    fourier transformation step.  The routine works in place. */
00100 void fourier::_fft_1d_2r (nr_double_t * r1, nr_double_t * r2, int len) {
00101   int n3, n2, j;
00102   nr_double_t rep, rem, aip, aim;
00103   n3 = 1 + (n2 = len + len);
00104 
00105   // put the two real vectors into one complex vector
00106   for (j = 1; j <= n2; j += 2) {
00107     r1[j] = r2[j-1];
00108   }
00109 
00110   // transform the complex vector
00111   _fft_1d (r1, len, 1);
00112 
00113   // separate the two transforms
00114   r2[0] = r1[1];
00115   r1[1] = r2[1] = 0.0;
00116   for (j = 2; j <= len; j += 2) {
00117     // use symmetries to separate the two transforms
00118     rep = 0.5 * (r1[j] + r1[n2-j]);
00119     rem = 0.5 * (r1[j] - r1[n2-j]);
00120     aip = 0.5 * (r1[j+1] + r1[n3-j]);
00121     aim = 0.5 * (r1[j+1] - r1[n3-j]);
00122     // ship them out in two complex vectors
00123     r1[j+1] = aim;
00124     r2[j+1] = -rem;
00125     r1[j] = r1[n2-j] = rep;
00126     r2[j] = r2[n2-j] = aip;
00127     r1[n3-j] = -aim;
00128     r2[n3-j] = rem;
00129   }
00130 }
00131 
00132 /* The following function transforms two vectors yielding real valued
00133    vectors using a single inverse fast fourier transformation step.
00134    The routine works in place as well. */
00135 void fourier::_ifft_1d_2r (nr_double_t * r1, nr_double_t * r2, int len) {
00136   nr_double_t r, i;
00137   int j, jj, nn = len + len;
00138 
00139   // put the two complex vectors into one complex vector
00140   for (j = 0, jj = 0; j < nn; j += 2) {
00141     r = r1[j] - r2[j+1];
00142     i = r1[j+1] + r2[j];
00143     r1[jj++] = r;
00144     r1[jj++] = i;
00145   }
00146 
00147   // transform the complex vector
00148   _fft_1d (r1, len, -1);
00149 
00150   // split the transforms into two real vectors
00151   for (j = 0; j < nn; j += 2) {
00152     r2[j] = r1[j+1];
00153     r1[j+1] = r2[j+1] = 0.0;
00154   }
00155 }
00156 
00157 /* This function performs a 1-dimensional fast fourier transformation
00158    on the given vector 'var'.  If 'sign' is -1 the inverse fft is
00159    computed, if +1 the fft itself is computed.  It returns a vector of
00160    binary size (as necessary for a fft). */
00161 vector fourier::fft_1d (vector var, int isign) {
00162   int i, n, len = var.getSize ();
00163 
00164   // compute necessary binary data array size
00165   int size = 2;
00166   while (size < len) size <<= 1;
00167 
00168   // put data items (temporary array) in place
00169   nr_double_t * data;
00170   data = (nr_double_t *) calloc (2 * size * sizeof (nr_double_t), 1);
00171   for (n = i = 0; i < len; i++, n += 2) {
00172     data[n] = real (var (i)); data[n+1] = imag (var (i));
00173   }
00174 
00175   // run 1-dimensional fft
00176   _fft_1d (data, size, isign);
00177 
00178   // store transformed data items in result vector
00179   vector res = vector (size);
00180   for (n = i = 0; i < size; i++, n += 2) {
00181     res (i) = nr_complex_t (data[n], data[n+1]);
00182     if (isign < 0) res (i) /= size;
00183   }
00184 
00185   // free temporary data array
00186   free (data);
00187   return res;
00188 }
00189 
00190 /* The function performs a 1-dimensional discrete fourier
00191    transformation.  Each data item is meant to be defined in
00192    equidistant steps. */
00193 void fourier::_dft_1d (nr_double_t * data, int len, int isign) {
00194   int k, n, size = 2 * len * sizeof (nr_double_t);
00195   nr_double_t * res = (nr_double_t *) calloc (size, 1);
00196   nr_double_t th, c, s;
00197   for (n = 0; n < 2 * len; n += 2) {
00198     th = n * pi / 2 / len;
00199     for (k = 0; k < 2 * len; k += 2) {
00200       c = cos (k * th);
00201       s = isign * sin (k * th);
00202       res[n] += data[k] * c + data[k+1] * s;
00203       res[n+1] += data[k+1] * c - data[k] * s;
00204     }
00205   }
00206   memcpy (data, res, size);
00207   free (res);
00208 }
00209 
00210 /* The function performs a 1-dimensional discrete fourier
00211    transformation on the given vector 'var'.  If 'sign' is -1 the
00212    inverse dft is computed, if +1 the dft itself is computed. */
00213 vector fourier::dft_1d (vector var, int isign) {
00214   int k, n, len = var.getSize ();
00215   vector res = vector (len);
00216   for (n = 0; n < len; n++) {
00217     nr_double_t th = - isign * 2 * pi * n / len;
00218     nr_complex_t val = 0;
00219     for (k = 0; k < len; k++)
00220       val += var (k) * std::polar (1.0, th * k);
00221     res (n) = isign < 0 ? val / (nr_double_t) len : val;
00222   }
00223   return res;
00224 }
00225 
00226 // Helper functions.
00227 vector fourier::ifft_1d (vector var) {
00228   return fft_1d (var, -1);
00229 }
00230 vector fourier::idft_1d (vector var) {
00231   return dft_1d (var, -1);
00232 }
00233 void fourier::_ifft_1d (nr_double_t * data, int len) {
00234   _fft_1d (data, len, -1);
00235 }
00236 void fourier::_idft_1d (nr_double_t * data, int len) {
00237   _dft_1d (data, len, -1);
00238 }
00239 
00240 /* The function performs a n-dimensional fast fourier transformation.
00241    Each data item is meant to be defined in equidistant steps.  The
00242    number of data items needs to be of binary size, e.g. 64, 128 for
00243    each dimension. */
00244 void fourier::_fft_nd (nr_double_t * data, int len[], int nd, int isign) {
00245 
00246  int i, i1, i2, i3, i2rev, i3rev, ip1, ip2, ip3, ifp1, ifp2;
00247  int ibit, k1, k2, n, np, nr, nt;
00248 
00249  // compute total number of complex values
00250  for (nt = 1, i = 0; i < nd; i++) nt *= len[i];
00251 
00252  // main loop over the dimensions
00253  for (np = 1, i = nd - 1; i >= 0; i--) {
00254    n = len[i];
00255    nr = nt / (n * np);
00256    ip1 = np << 1;
00257    ip2 = ip1 * n;
00258    ip3 = ip2 * nr;
00259 
00260    // bit-reversal method
00261    for (i2rev = 1, i2 = 1; i2 <= ip2; i2 += ip1) {
00262      if (i2 < i2rev) {
00263        for (i1 = i2; i1 <= i2 + ip1 - 2; i1 += 2) {
00264          for (i3 = i1; i3 <= ip3; i3 += ip2) {
00265            i3rev = i2rev + i3 - i2;
00266            Swap (data[i3-1], data[i3rev-1]);
00267            Swap (data[i3], data[i3rev]);
00268          }
00269        }
00270      }
00271      ibit = ip2 >> 1;
00272      while (ibit >= ip1 && i2rev > ibit) {
00273        i2rev -= ibit;
00274        ibit >>= 1;
00275      }
00276      i2rev += ibit;
00277    }
00278 
00279    // Danielson-Lanzcos algorithm
00280    nr_double_t ti, tr, wt, th, wr, wi, wpi, wpr;
00281    ifp1 = ip1;
00282    while (ifp1 < ip2) {
00283      ifp2 = ifp1 << 1;
00284      th = isign * 2 * pi / (ifp2 / ip1);
00285      wt = sin (0.5 * th);
00286      wpr = -2.0 * wt * wt;
00287      wpi = sin (th);
00288      wr = 1.0;
00289      wi = 0.0;
00290      for (i3 = 1; i3 <= ifp1; i3 += ip1) {
00291        for (i1 = i3; i1 <= i3 + ip1 - 2; i1 += 2) {
00292          for (i2 = i1; i2 <= ip3; i2 += ifp2) {
00293            k1 = i2;
00294            k2 = k1 + ifp1;
00295            tr = wr * data[k2-1] - wi * data[k2];
00296            ti = wr * data[k2] + wi * data[k2-1];
00297            data[k2-1] = data[k1-1] - tr;
00298            data[k2] = data[k1] - ti;
00299            data[k1-1] += tr;
00300            data[k1] += ti;
00301          }
00302        }
00303        wr = (wt = wr) * wpr - wi * wpi + wr;
00304        wi = wi * wpr + wt * wpi + wi;
00305      }
00306      ifp1 = ifp2;
00307    }
00308    np *= n;
00309  }
00310 }
00311 
00312 // Helper functions.
00313 void fourier::_ifft_nd (nr_double_t * data, int len[], int nd) {
00314   _fft_nd (data, len, nd, -1);
00315 }
00316 
00317 // Shuffles values of FFT around.
00318 vector fourier::fftshift (vector var) {
00319   int i, n, len = var.getSize ();
00320   vector res = vector (len);
00321   n = len / 2;
00322   for (i = 0; i < len / 2; i++) {
00323     res (i) = var (n + i);
00324     res (i + n) = var (i);
00325   }
00326   return res;
00327 }
00328 
00329 } // namespace qucs