htswanalysis/MACS/lib/gsl/gsl-1.11/fft/signals_source.c

   1 /* fft/signals_source.c
   2  *
   3  * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2007 Brian Gough
   4  *
   5  * This program is free software; you can redistribute it and/or modify
   6  * it under the terms of the GNU General Public License as published by
   7  * the Free Software Foundation; either version 3 of the License, or (at
   8  * your option) any later version.
   9  *
  10  * This program is distributed in the hope that it will be useful, but
  11  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13  * General Public License for more details.
  14  *
  15  * You should have received a copy of the GNU General Public License
  16  * along with this program; if not, write to the Free Software
  17  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
  18  */
  19
  20 #include "signals.h"
  21
  22 int
  23 FUNCTION(fft_signal,complex_pulse) (const size_t k,
  24                                     const size_t n,
  25                                     const size_t stride,
  26                                     const BASE z_real,
  27                                     const BASE z_imag,
  28                                     BASE data[],
  29                                     BASE fft[])
  30 {
  31   size_t j;
  32
  33   if (n == 0)
  34     {
  35       GSL_ERROR ("length n must be positive integer", GSL_EDOM);
  36     }
  37
  38   /* gsl_complex pulse at position k,  data[j] = z * delta_{jk} */
  39
  40   for (j = 0; j < n; j++)
  41     {
  42       REAL(data,stride,j) = 0.0;
  43       IMAG(data,stride,j) = 0.0;
  44     }
  45
  46   REAL(data,stride,k % n) = z_real;
  47   IMAG(data,stride,k % n) = z_imag;
  48
  49   /* fourier transform, fft[j] = z * exp(-2 pi i j k / n) */
  50
  51   for (j = 0; j < n; j++)
  52     {
  53       const double arg = -2 * M_PI * ((double) ((j * k) % n)) / ((double) n);
  54       const BASE w_real = (BASE)cos (arg);
  55       const BASE w_imag = (BASE)sin (arg);
  56       REAL(fft,stride,j) = w_real * z_real - w_imag * z_imag;
  57       IMAG(fft,stride,j) = w_real * z_imag + w_imag * z_real;
  58     }
  59
  60   return 0;
  61
  62 }
  63
  64
  65 int
  66 FUNCTION(fft_signal,complex_constant) (const size_t n,
  67                                        const size_t stride,
  68                                        const BASE z_real,
  69                                        const BASE z_imag,
  70                                        BASE data[],
  71                                        BASE fft[])
  72 {
  73   size_t j;
  74
  75   if (n == 0)
  76     {
  77       GSL_ERROR ("length n must be positive integer", GSL_EDOM);
  78     }
  79
  80   /* constant, data[j] = z */
  81
  82   for (j = 0; j < n; j++)
  83     {
  84       REAL(data,stride,j) = z_real;
  85       IMAG(data,stride,j) = z_imag;
  86     }
  87
  88   /* fourier transform, fft[j] = n z delta_{j0} */
  89
  90   for (j = 0; j < n; j++)
  91     {
  92       REAL(fft,stride,j) = 0.0;
  93       IMAG(fft,stride,j) = 0.0;
  94     }
  95
  96   REAL(fft,stride,0) = ((BASE) n) * z_real;
  97   IMAG(fft,stride,0) = ((BASE) n) * z_imag;
  98
  99   return 0;
 100
 101 }
 102
 103
 104 int
 105 FUNCTION(fft_signal,complex_exp) (const int k,
 106                                   const size_t n,
 107                                   const size_t stride,
 108                                   const BASE z_real,
 109                                   const BASE z_imag,
 110                                   BASE data[],
 111                                   BASE fft[])
 112 {
 113   size_t j;
 114
 115   if (n == 0)
 116     {
 117       GSL_ERROR ("length n must be positive integer", GSL_EDOM);
 118     }
 119
 120   /* exponential,  data[j] = z * exp(2 pi i j k) */
 121
 122   for (j = 0; j < n; j++)
 123     {
 124       const double arg = 2 * M_PI * ((double) ((j * k) % n)) / ((double) n);
 125       const BASE w_real = (BASE)cos (arg);
 126       const BASE w_imag = (BASE)sin (arg);
 127       REAL(data,stride,j) = w_real * z_real - w_imag * z_imag;
 128       IMAG(data,stride,j) = w_real * z_imag + w_imag * z_real;
 129     }
 130
 131   /* fourier transform, fft[j] = z * delta{(j - k),0} */
 132
 133   for (j = 0; j < n; j++)
 134     {
 135       REAL(fft,stride,j) = 0.0;
 136       IMAG(fft,stride,j) = 0.0;
 137     }
 138
 139   {
 140     int freq;
 141
 142     if (k <= 0)
 143       {
 144         freq = (n-k) % n ;
 145       }
 146     else
 147       {
 148         freq = (k % n);
 149       };
 150
 151     REAL(fft,stride,freq) = ((BASE) n) * z_real;
 152     IMAG(fft,stride,freq) = ((BASE) n) * z_imag;
 153   }
 154
 155   return 0;
 156
 157 }
 158
 159
 160 int
 161 FUNCTION(fft_signal,complex_exppair) (const int k1,
 162                                       const int k2,
 163                                       const size_t n,
 164                                       const size_t stride,
 165                                       const BASE z1_real,
 166                                       const BASE z1_imag,
 167                                       const BASE z2_real,
 168                                       const BASE z2_imag,
 169                                       BASE data[],
 170                                       BASE fft[])
 171 {
 172   size_t j;
 173
 174   if (n == 0)
 175     {
 176       GSL_ERROR ("length n must be positive integer", GSL_EDOM);
 177     }
 178
 179   /* exponential,  data[j] = z1 * exp(2 pi i j k1) + z2 * exp(2 pi i j k2) */
 180
 181   for (j = 0; j < n; j++)
 182     {
 183       const double arg1 = 2 * M_PI * ((double) ((j * k1) % n)) / ((double) n);
 184       const BASE w1_real = (BASE)cos (arg1);
 185       const BASE w1_imag = (BASE)sin (arg1);
 186       const double arg2 = 2 * M_PI * ((double) ((j * k2) % n)) / ((double) n);
 187       const BASE w2_real = (BASE)cos (arg2);
 188       const BASE w2_imag = (BASE)sin (arg2);
 189       REAL(data,stride,j) = w1_real * z1_real - w1_imag * z1_imag;
 190       IMAG(data,stride,j) = w1_real * z1_imag + w1_imag * z1_real;
 191       REAL(data,stride,j) += w2_real * z2_real - w2_imag * z2_imag;
 192       IMAG(data,stride,j) += w2_real * z2_imag + w2_imag * z2_real;
 193     }
 194
 195   /* fourier transform, fft[j] = z1 * delta{(j - k1),0} + z2 *
 196      delta{(j - k2,0)} */
 197
 198   for (j = 0; j < n; j++)
 199     {
 200       REAL(fft,stride,j) = 0.0;
 201       IMAG(fft,stride,j) = 0.0;
 202     }
 203
 204   {
 205     int freq1, freq2;
 206
 207     if (k1 <= 0)
 208       {
 209         freq1 = (n - k1) % n;
 210       }
 211     else
 212       {
 213         freq1 = (k1 % n);
 214       };
 215
 216     if (k2 <= 0)
 217       {
 218         freq2 = (n - k2) % n;
 219       }
 220     else
 221       {
 222         freq2 = (k2 % n);
 223       };
 224
 225     REAL(fft,stride,freq1) += ((BASE) n) * z1_real;
 226     IMAG(fft,stride,freq1) += ((BASE) n) * z1_imag;
 227     REAL(fft,stride,freq2) += ((BASE) n) * z2_real;
 228     IMAG(fft,stride,freq2) += ((BASE) n) * z2_imag;
 229   }
 230
 231   return 0;
 232
 233 }
 234
 235
 236 int
 237 FUNCTION(fft_signal,complex_noise) (const size_t n,
 238                                     const size_t stride,
 239                                     BASE data[],
 240                                     BASE fft[])
 241 {
 242   size_t i;
 243   int status;
 244
 245   if (n == 0)
 246     {
 247       GSL_ERROR ("length n must be positive integer", GSL_EDOM);
 248     }
 249
 250   for (i = 0; i < n; i++)
 251     {
 252       REAL(data,stride,i) = (BASE)urand();
 253       IMAG(data,stride,i) = (BASE)urand();
 254     }
 255
 256   /* compute the dft */
 257   status = FUNCTION(gsl_dft_complex,forward) (data, stride, n, fft);
 258
 259   return status;
 260 }
 261
 262
 263 int
 264 FUNCTION(fft_signal,real_noise) (const size_t n,
 265                                  const size_t stride,
 266                                  BASE data[],
 267                                  BASE fft[])
 268 {
 269   size_t i;
 270   int status;
 271
 272   if (n == 0)
 273     {
 274       GSL_ERROR ("length n must be positive integer", GSL_EDOM);
 275     }
 276
 277   for (i = 0; i < n; i++)
 278     {
 279       REAL(data,stride,i) = (BASE)urand();
 280       IMAG(data,stride,i) = 0.0;
 281     }
 282
 283   /* compute the dft */
 284   status = FUNCTION(gsl_dft_complex,forward) (data, stride, n, fft);
 285
 286   return status;
 287 }
 288