1 /* fft/signals_source.c
3 * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2007 Brian Gough
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.
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.
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.
23 FUNCTION(fft_signal,complex_pulse) (const size_t k,
35 GSL_ERROR ("length n must be positive integer", GSL_EDOM);
38 /* gsl_complex pulse at position k, data[j] = z * delta_{jk} */
40 for (j = 0; j < n; j++)
42 REAL(data,stride,j) = 0.0;
43 IMAG(data,stride,j) = 0.0;
46 REAL(data,stride,k % n) = z_real;
47 IMAG(data,stride,k % n) = z_imag;
49 /* fourier transform, fft[j] = z * exp(-2 pi i j k / n) */
51 for (j = 0; j < n; j++)
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;
66 FUNCTION(fft_signal,complex_constant) (const size_t n,
77 GSL_ERROR ("length n must be positive integer", GSL_EDOM);
80 /* constant, data[j] = z */
82 for (j = 0; j < n; j++)
84 REAL(data,stride,j) = z_real;
85 IMAG(data,stride,j) = z_imag;
88 /* fourier transform, fft[j] = n z delta_{j0} */
90 for (j = 0; j < n; j++)
92 REAL(fft,stride,j) = 0.0;
93 IMAG(fft,stride,j) = 0.0;
96 REAL(fft,stride,0) = ((BASE) n) * z_real;
97 IMAG(fft,stride,0) = ((BASE) n) * z_imag;
105 FUNCTION(fft_signal,complex_exp) (const int k,
117 GSL_ERROR ("length n must be positive integer", GSL_EDOM);
120 /* exponential, data[j] = z * exp(2 pi i j k) */
122 for (j = 0; j < n; j++)
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;
131 /* fourier transform, fft[j] = z * delta{(j - k),0} */
133 for (j = 0; j < n; j++)
135 REAL(fft,stride,j) = 0.0;
136 IMAG(fft,stride,j) = 0.0;
151 REAL(fft,stride,freq) = ((BASE) n) * z_real;
152 IMAG(fft,stride,freq) = ((BASE) n) * z_imag;
161 FUNCTION(fft_signal,complex_exppair) (const int k1,
176 GSL_ERROR ("length n must be positive integer", GSL_EDOM);
179 /* exponential, data[j] = z1 * exp(2 pi i j k1) + z2 * exp(2 pi i j k2) */
181 for (j = 0; j < n; j++)
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;
195 /* fourier transform, fft[j] = z1 * delta{(j - k1),0} + z2 *
198 for (j = 0; j < n; j++)
200 REAL(fft,stride,j) = 0.0;
201 IMAG(fft,stride,j) = 0.0;
209 freq1 = (n - k1) % n;
218 freq2 = (n - k2) % n;
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;
237 FUNCTION(fft_signal,complex_noise) (const size_t n,
247 GSL_ERROR ("length n must be positive integer", GSL_EDOM);
250 for (i = 0; i < n; i++)
252 REAL(data,stride,i) = (BASE)urand();
253 IMAG(data,stride,i) = (BASE)urand();
256 /* compute the dft */
257 status = FUNCTION(gsl_dft_complex,forward) (data, stride, n, fft);
264 FUNCTION(fft_signal,real_noise) (const size_t n,
274 GSL_ERROR ("length n must be positive integer", GSL_EDOM);
277 for (i = 0; i < n; i++)
279 REAL(data,stride,i) = (BASE)urand();
280 IMAG(data,stride,i) = 0.0;
283 /* compute the dft */
284 status = FUNCTION(gsl_dft_complex,forward) (data, stride, n, fft);