#define PROG_NAME "blimpulse" #define PROG_DESC "compute band-limited finite pulses" #define PROG_VERS "1.1" /* Copyright © 2007 by the State University of Campinas (UNICAMP). */ /* Last edited on 2024-12-01 00:26:40 by stolfi */ #define _GNU_SOURCE #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include /* GENERAL PARAMETERS */ /* INTERNAL PROTOTYPES */ int main (int argc, char **argv); void find_compact_pulse(int ns, int nt, int type, double X[]); /* Tries to find a pulse with {nt} samples that is compact in time (many consecutive zeros) and frequency (many highest freq components are absent). The nonzero samples of the pulse will be contained in {X[0..ns-1]}, which are set by the procedure. Samples {X[ns..nt-1]} are implicitly zero. The {type} argument chooses between various normalization constraints. */ void fourier_to_spectrum(int n, fftw_complex ot[], double P[]); /* Reduces a complex Fourier transform {ot[0..n-1][0..1]} into a power spectrum {P[0..fMax]} where {fMax == floor(n/2)}. Also normalizes the result so that the total energy (Euclidean norm squared) of the original signal is the sum of all entries in {P}. */ void pulse_spectrum ( int ns, int nt, double X[], fftw_complex in[], fftw_complex ot[], fftw_plan plan, double fW[] ); /* Pads the pulse {X[0..ns-1]} with zeros to width {nt}, then computes its FFT and reduces it to the extended power spectrum {P[0..fMax]} where {fMax == nt/2}. Assumes that the vectors {in} and {ot} have {nt} complex elements. */ double pulse_badness ( int ns, int nt, double X[], fftw_complex in[], fftw_complex ot[], fftw_plan plan, double fW[] ); /* A badness score that measures the high freq contents of the pulse {X[0..ns-1]} after padding to {nt} samples. Assumes that {in} and {ot} have {nt} complex elements. */ /* IMPLEMENTATIONS */ int main (int argc, char **argv) { demand (argc == 4, "usage: %s {NTOTAL} {NNONZ} {NORMTYPE}"); int nt = atoi(argv[1]); /* Number of samples in period. */ int ns = atoi(argv[2]); /* Number of possibly nonzero samples in pulse. */ int type = atoi(argv[3]); /* Normalization type. */ double X[ns]; find_compact_pulse(ns, nt, type, X); fclose(stderr); fclose(stdout); return (0); } void find_compact_pulse(int ns, int nt, int type, double X[]) { bool_t debug = FALSE; char *fprefix = jsprintf("out/ot-%04d-%04d-%d", nt, ns, type); /* Working storage for the goal function: */ fftw_complex *in = fftw_malloc(sizeof(fftw_complex)*nt); fftw_complex *ot = fftw_malloc(sizeof(fftw_complex)*nt); fftw_plan plan = fftw_plan_dft_1d(nt, in, ot, FFTW_FORWARD, FFTW_ESTIMATE); /* Compute the max freq for {nt} samples: */ demand((nt % 2) == 1, "{nt} must be odd"); int fMax = (nt-1)/2; /* Compute the max freq allowed {fOkMax}: */ demand(ns < nt, "{ns} must be less than {nt}"); demand((ns % 2) == 1, "{ns} must be odd"); int fOkMax = (nt - ns)/2; assert(fOkMax > 0); /* Compute the number {nk} of freqs to kill: */ int nk = ns - 1; /* We should have just one freedom left. */ /* For example if {nt = 17} and {ns = 7}, then {fOkMax} will be {(17 - 7)/2 = 5} and {nk} will be 6, meaning that the 6 frequencies {-8,-7,-6,+6,+7,+8} outside the range {-5..+5} will be suppressed. */ /* Build the linear system that says that hi-freq terms are zero. */ int rows = nk + 1; /* {nk} freq-kill eqs plus one unit-mag eq. */ int cols = ns; /* {ns} unknowns (the samples {X[0..ns-1]}). */ assert(cols == rows); /* The system should be square. */ double A[rows*cols]; /* The constraint coefficient matrix. */ double B[rows]; /* The independent-term vector. */ double R[rows]; /* Residuals of constraints. */ double Y[cols]; /* Estimated error in {X} computed from residuals. */ /* The {nk} frequency-kill constraints: */ int f; int ie = 0; for (f = fOkMax+1; f <= fMax; f++) { int sgn; for (sgn = -1; sgn <= +1; sgn += 2) { /* Require that frequency {f*sgn} be absent: */ int j; for (j = 0; j < ns; j++) { A[ie*cols + j] = sin(M_PI/4 + (2*M_PI*f*sgn*j)/nt); } B[ie] = 0; ie++; } } assert(ie == nk); /* Unit-sum constraint: */ int iu = nk; /* Row of the unit-mag constraint in {A} */ if (type == 0) { /* Require the center sample to be 1: */ int j; for (j = 0; j < ns; j++) { A[iu*cols + j] = 0; } A[iu*cols + ns/2] = 1; } else if (type == 1) { /* Require slope 1 at the center sample. */ int j; for (j = 0; j < ns; j++) { A[iu*cols + j] = 0; } A[iu*cols + ns/2-1] = -0.5; A[iu*cols + ns/2+1] = +0.5; } else if (type == 2) { /* Require the sum of all samples to be 1. */ int j; for (j = 0; j < ns; j++) { A[iu*cols + j] = 1; } } else { demand(FALSE, "invalid normalization {type}"); } B[iu] = 1; /* Solve the system: */ uint32_t r1; gausol_solve(rows, cols, A, 1, B, X, TRUE,TRUE, 0.0, NULL, &r1); demand(r1 == rows, "indeterminate system"); int iter; for (iter = 0; iter < 20; iter++) { /* Now {X} is a tentative solution: */ if (debug) { rn_gen_print(stderr, ns, X, "%20.16f", "X =\n ", "\n ", "\n"); } /* Compute the residual error {R}: */ int r2 = gausol_residual(rows, cols, A, 1, B, X, R); demand(r2 == rows, "indeterminate system"); if (debug) { rn_gen_print(stderr, rows, R, "%20.16f", "R =\n ", "\n ", "\n"); } /* Solve the residual system to get a correction {Y}: */ uint32_t r3; gausol_solve(rows, cols, A, 1, R, Y, TRUE,TRUE, 0.0, NULL, &r3); demand(r3 == rows, "indeterminate system"); if (debug) { rn_gen_print(stderr, ns, Y, "%20.16f", "Y =\n ", "\n ", "\n"); } double Rnorm = rn_L_inf_norm(rows, R); double Ynorm = rn_L_inf_norm(cols, Y); fprintf(stderr, " iteration %4d residual norm = %20.16f correction norm = %20.16f\n", iter, Rnorm, Ynorm); /* Subtract correction from solution: */ rn_sub(cols, X, Y, X); } /* Print final solution: */ rn_gen_print(stderr, ns, X, "%20.16f", "X =\n ", "\n ", "\n"); /* Write the pulse for plotting: */ { FILE *wr = open_write(txtcat(fprefix,"-pulse.dat"), TRUE); int i; double Lipp = INFINITY; /* {log(X[i-2])} */ double Lip = INFINITY; /* {log(X[i-1])} */ for (i = -2; i < nt+2; i++) { int k = ((i % nt) + nt) % nt; double Xi = (k < ns ? X[k] : 0); double Yi = (k < ns ? Y[k] : 0); /* Logs and derivatives of logs: */ double Li = log((fabs(Xi) <= 0 ? 1e-300 : fabs(Xi))); double DLi = (Lip == INFINITY ? 0 : Li - Lip); double DDLi = (Lipp == INFINITY ? 0 : (Li - 2*Lip + Lipp)/2); fprintf ( wr, "%5d %20.16f %+12.5f %+12.5f %+12.5f %20.16f\n", i, Xi, Li, DLi, DDLi, Yi ); Lipp = Lip; Lip = Li; } fclose(wr); } /* Write the power spectrum for plotting: */ { double P[fMax+1]; pulse_spectrum(ns, nt, X, in, ot, plan, P); FILE *wr = open_write(txtcat(fprefix,"-spectrum.dat"), TRUE); int f; double Lfpp = INFINITY; /* {log(P[f-2])} */ double Lfp = INFINITY; /* {log(P[f-1])} */ for (f = -fMax; f <= +fMax; f++) { /* Use the bilateral spectrum, it plots nicer than the folded one: */ double Pf = (f == 0 ? P[f] : P[abs(f)]/2); /* Use the bilateral spectrum, it plots nicer than the folded one: */ /* Logs and derivatives of logs: */ double Lf = log((Pf <= 0 ? 1e-300 : Pf)); double DLf = (Lfp == INFINITY ? 0 : Lf - Lfp); double DDLf = (Lfpp == INFINITY ? 0 : (Lf - 2*Lfp + Lfpp)/2); fprintf ( wr, "%5d %20.16f %+12.5f %+12.5f %+12.5f\n", f, Pf, Lf, DLf, DDLf ); Lfpp = Lfp; Lfp = Lf; } fclose(wr); } } void pulse_spectrum ( int ns, int nt, double X[], fftw_complex in[], fftw_complex ot[], fftw_plan plan, double P[] ) { /* Compute the Fourier transform: */ int i; for (i = 0; i < nt; i++) { in[i][0] = (i < ns ? X[i] : 0); in[i][1] = 0; } fftw_execute(plan); /* Extract the power spectrum from FFT. */ fourier_to_spectrum(nt, ot, P); } void fourier_to_spectrum(int nt, fftw_complex ot[], double P[]) { int fMax = nt/2; /* Arg shoud be {fftw_complex} not {double *} but gcc complains. */ auto double cmod2(double *c); /* Complex modulus of {c}, squared. */ double cmod2(double *c) { return c[0]*c[0] + c[1]*c[1]; } /* We must combine Fourier terms with same freq. */ /* They are {ot[i]} and {ot[n-i]}, modulo {n}. */ /* Frequency 0 has a single term: */ P[0] = cmod2(ot[0])/nt; /* Other freqs have two terms, except for {fMax} when {nt} is even: */ int f; for (f = 1; f <= fMax; f++) { double Pi = cmod2(ot[f]); int j = nt - f; if (j != f) { Pi += cmod2(ot[j]); } P[f] = Pi/nt; } }