Qucs-core
0.0.19
|
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