Code:
#include <stdlib.h>
#include <stdint.h>
#include <string.h>
#include <stdio.h>
#include <math.h>
#include <complex.h>
#include <fftw3.h>
static uint64_t prng64(uint64_t x)
{
x ^= x >> 12;
x ^= x << 25;
x ^= x >> 27;
return x * (uint64_t)2685821657736338717ULL;
}
static double f_noise(const double v)
{
uint64_t x = (uint64_t)(1099511627776.0 * v);
int n = 16;
while (n-->0)
x = prng64(x);
return 1.0 - (x & (uint64_t)2199023255551ULL) / 1099511627776.0;
}
static double f_default(const double x)
{
const double r = x * 3.14159265358979323864;
return sin(r) * sin(20.0*r);
}
static double f_sawpulse(const double x)
{
const double w = 16.0 * x;
return x*x*(x*x - x - x + 1.0) * 32.0 * ((w - (long)w) - 0.5);
}
static double f_pulse(const double x)
{
const double r = x * 3.14159265358979323864;
return sin(r)*sin(r) * sin(15.1*r);
}
int main(int argc, char *argv[])
{
double (*sample)(const double) = f_default;
double *s1, *sn;
fftw_complex *c1, *cn;
fftw_plan p;
long n, w, c, i;
char dummy;
if (argc < 3 || argc > 4) {
fprintf(stderr, "\n");
fprintf(stderr, "Usage: %s SAMPLES WAVEPACKETS [ WAVEFORM ]\n", argv[0]);
fprintf(stderr, "\n");
fprintf(stderr, "Waveforms: (t=0..1)\n");
fprintf(stderr, " default sin(pi*t)*sin(20*pi*t)\n");
fprintf(stderr, " pulse sin(2*pi*t)*sin(2*pi*t)*sin(15.1*pi*t)\n");
fprintf(stderr, " sawpulse t*t*(16*t*t-32*t+16)*((16*t) mod 1.0)\n");
fprintf(stderr, " noise Xorshift*-based noise\n");
fprintf(stderr, "\n");
return EXIT_FAILURE;
}
if (argc > 3) {
if (!strcmp(argv[3], "sawpulse"))
sample = f_sawpulse;
else
if (!strcmp(argv[3], "noise"))
sample = f_noise;
else
if (!strcmp(argv[3], "pulse"))
sample = f_pulse;
else
if (!strcmp(argv[3], "default"))
sample = f_default;
else {
fprintf(stderr, "%s: Unknown waveform.\n", argv[3]);
return EXIT_FAILURE;
}
}
if (sscanf(argv[1], " %ld %c", &n, &dummy) != 1 || n < 1L) {
fprintf(stderr, "%s: Invalid number of samples.\n", argv[1]);
return EXIT_FAILURE;
}
if (sscanf(argv[2], " %ld %c", &w, &dummy) != 1 || w < 1L || w >= n) {
fprintf(stderr, "%s: Invalid number of wave packets.\n", argv[2]);
return EXIT_FAILURE;
}
if (n % w) {
fprintf(stderr, "The number of samples is not a multiple of wave packets.\n");
return EXIT_FAILURE;
} else
c = n / w;
s1 = fftw_malloc((size_t)c * sizeof *s1);
sn = fftw_malloc((size_t)n * sizeof *sn);
c1 = fftw_malloc((size_t)c * sizeof *c1);
cn = fftw_malloc((size_t)n * sizeof *cn);
if (!s1 || !sn || !c1 || !cn) {
fprintf(stderr, "Not enough memory.\n");
return EXIT_FAILURE;
}
/* Initialize one wave packet to s1 and sn. */
for (i = 0L; i < c; i++)
s1[i] = sn[i] = sample((double)i / (double)c);
/* sn contains zeroes following the first wave packet. */
for (i = c; i < n; i++)
sn[i] = 0.0;
/* Clear c1 and cn. */
for (i = 0; i < c; i++)
c1[i] = 0.0;
for (i = 0; i < n; i++)
cn[i] = 0.0;
p = fftw_plan_dft_r2c_1d(c, s1, c1, FFTW_ESTIMATE);
fftw_execute(p);
fftw_destroy_plan(p);
p = fftw_plan_dft_r2c_1d(n, sn, cn, FFTW_ESTIMATE);
fftw_execute(p);
fftw_destroy_plan(p);
/* Fill in c1 and cn. */
for (i = 1L; i < c-i; i++)
c1[c - i] = conj(c1[i]);
for (i = 1L; i < n-i; i++)
cn[n - i] = conj(cn[i]);
printf("#i f[i] |F%ld[i]| Re(F%ld[i]) Im(F%ld[i]) |F[i]| Re(F[i]) Im(F[i]) i/%ld\n", w, w, w, w);
for (i = 0; i < n; i++)
if (i % w == 0)
printf("%.9f %.9f %.9f %.9f %.9f %.9f %.9f %.9f %.9f\n",
(double)i / (double)n, sn[i],
cabs(cn[i]), creal(cn[i]), cimag(cn[i]),
cabs(c1[i/w]), creal(c1[i/w]), cimag(c1[i/w]), (double)(i/w) / (double)n);
else
printf("%.9f %.9f %.9f %.9f %.9f\n",
(double)i / (double)n, sn[i],
cabs(cn[i]), creal(cn[i]), cimag(cn[i]));
return EXIT_SUCCESS;
}
The