/* Last edited on 2005-01-01 11:31:53 by stolfi */ double frb_mean_dist_sqr (frb_curve_t *a, frb_curve_t *b); /* Returns the average distance (squared) between coresponding samples of the two curves. */ double frb_mean_dist_sqr (frb_curve_t *a, frb_curve_t *b) { double d2 = 0.0; int NA = a->nel; int NB = b->nel; int ns = (NA > NB ? NA : NB); d2 = r3_dist_sqr(&(a->el[0]), &(b->el[0])); int i; for (i = 1; i < ns; i++) { int i0 = (int)(0.5 + ((double)(NA-1)*i)/((double)(ns-1))); int i1 = (int)(0.5 + ((double)(NB-1)*i)/((double)(ns-1))); d2 += r3_dist_sqr(&(a->el[i0]), &(b->el[i1])); } return 2.0 * d2/((double)(NA + NB)); } i2_t frb_find_best_subcand_quick ( frb_curve_t *a, frb_curve_t *b, i2_t ini, i2_t ctr, i2_t fin, int nSteps, int maxAdjust ); /* Given two curves {a,b}, finds the pair of indices {ia,ib} such that matching the sub-segment {a[ia-h..ia+h]} with {b[ib-h..ib+h} minimizes the average square distance between samples, where {h = nSteps/2}. Considers only {2*maxAdjust+1} possible segments. After pairing the {a} and {b} samples by the trivial (``most perfect'') matching, considers all choices of {ia} and {ib} that are symmetrically displaced from the centers {aCtr,bCtr} of those segments, and whose magnitude is at most {maxAdjust}. */ i2_t frb_find_best_subcand_quick ( frb_curve_t *a, frb_curve_t *b, i2_t ini, i2_t ctr, i2_t fin, int nSteps, int maxAdjust ) { affirm(nSteps % 2 == 0 , "nSteps must be even"); i2_t ctrBest; double d2Best = INF; int shiftBest = 0; int hSteps = nSteps / 2; if (debug_align) { fprintf(stderr, "# aligning a[%d..%d] b[%d..%d]", ini.c[0], fin.c[0], ini.c[1], fin.c[1] ); } int shift; for (shift = -maxAdjust; shift <= maxAdjust; shift++) { int da = shift / 2; int db = shift - da; int aCtr = ctr.c[0] + da, aIni = aCtr - hSteps, aFin = aCtr + hSteps; int bCtr = ctr.c[1] + db, bIni = bCtr - hSteps, bFin = bCtr + hSteps; if ( (aIni >= ini.c[0]) && (aFin <= fin.c[0]) && (bIni >= ini.c[1]) && (bFin <= fin.c[1]) ) { if (debug_align) { fprintf(stderr, "# shift = %3d a = [%d..%d] b = [%d..%d]", shift, aIni, aFin, bIni, bFin ); } frb_curve_t ca = frb_curve_trim(a, aIni, nSteps+1); frb_curve_t cb = frb_curve_trim(b, bIni, nSteps+1); r4x4_t m = frb_curve_alignment_matrix(&ca, &cb, 0.20); frb_curve_t cam = frb_curve_map(&ca, &m); double d2 = frb_mean_dist_sqr(&cam, &cb); if (debug_align) { fprintf(stderr, " dist = %8.3f\n", sqrt(d2)); } if (d2 < d2Best) { d2Best = d2; shiftBest = shift; ctrBest = (i2_t){{ aCtr, bCtr }}; } free(cam.el); } } fprintf(stderr, "# best shift = %3d", shiftBest); fprintf(stderr, " dist = %8.3f (best)\n", sqrt(d2Best)); if (debug_align) { affirm(FALSE, "debug on - aborted"); } return ctrBest; } i2_t frb_find_best_subcand_slow ( frb_curve_t *a, frb_curve_t *b, i2_t ini, i2_t ctr, i2_t fin, int nSteps, int maxAdjust ); /* Same as {frb_find_best_subcand_quick}, but considers all sub-segments of {a} and {b} that are displaced by at most {maxShift} from the trivial alignment. */ i2_t frb_find_best_subcand_slow ( frb_curve_t *a, frb_curve_t *b, i2_t ini, i2_t ctr, i2_t fin, int nSteps, int maxAdjust ) { affirm(nSteps % 2 == 0 , "nSteps must be even"); int shiftBest = 0; i2_t ctrBest; double d2Best = INF; int hSteps = nSteps / 2; if (debug_align) { fprintf(stderr, "# aligning a[%d..%d] b[%d..%d]", ini.c[0], fin.c[0], ini.c[1], fin.c[1] ); } int shift; for (shift = -maxAdjust; shift <= maxAdjust; shift++) { int da = shift / 2, aMid = ctr.c[0] - da; int db = shift - da, bMid = ctr.c[1] + db; int aLo = ini.c[0] + hSteps - aMid, aHi = fin.c[0] - hSteps - aMid; int bLo = ini.c[1] + hSteps - bMid, bHi = fin.c[1] - hSteps - bMid; int mLo = (aLo > bLo ? aLo : bLo), mHi = (aHi < bHi ? aHi : bHi); if (debug_align) { fprintf(stderr, "# shift = %d aMid = %3d bMid = %3d", shift, aMid, bMid); fprintf(stderr, " aLo = %3d aHi = %3d", aLo, aHi); fprintf(stderr, " bLo = %3d bHi = %3d", bLo, bHi); fprintf(stderr, " mLo = %3d mHi = %3d", mLo, mHi); fprintf(stderr, "\n"); } int mc; for (mc = mLo; mc <= mHi; mc++) { int aCtr = aMid + mc; int bCtr = bMid + mc; if (debug_align) { fprintf(stderr, "# aCtr = %3d bCtr = %3d", aCtr, bCtr); } double d2 = frb_segment_discepancy(a, aCtr, b bCtr, nSteps); if (debug_align) { fprintf(stderr, " dist = %8.3f\n", sqrt(d2)); } if (d2 < d2Best) { d2Best = d2; shiftBest = shift; ctrBest = (i2_t){{ aCtr, bCtr }}; } } } fprintf(stderr, "# best shift = %3d aCtr = %d bCtr = %d", shiftBest, ctrBest.c[0], ctrBest.c[1] ); fprintf(stderr, " dist = %8.3f (best)\n", sqrt(d2Best)); if (debug_align) { affirm(FALSE, "debug on - aborted"); } return ctrBest; } frb_curve_t frb_pack_stats (frb_transform_t *avg, frb_transform_t *var); /* Given the sample-wise averages {avg[i]} and variances {var[i]} of a set of transforms, returns a vector of {r3_t}s such that {r[i].c[0..2]} are the average, lower bound, and upper bound of the corresponding sample, where the bounds are the average plus or minus one standard deviation. */ frb_curve_t frb_pack_stats (frb_transform_t *avg, frb_transform_t *var) { int fmax = avg->nel; affirm(fmax == var->nel, "inconsistent lengths"); frb_curve_t res = frb_curve_new(fmax); int p; for (p = 0; p < fmax; p++) { double ap = avg->el[p]; double vp = var->el[p]; double dp = sqrt(vp); r3_t *rp = &(res.el[p]); rp->c[0] = ap; rp->c[1] = ap - dp; rp->c[2] = ap + dp; } return res; } bool_t logScale; /* TRUE computes the log-mean-exp of the signals. */ o->logScale = argparser_keyword_present(pp, "-logScale"); void frb_to_log (double eps, frb_signal_t *s); /* Replaces each element {s[i]} by the log of its absolute value, biased by the noise-level {eps} -- that is, {0.5*log(eps^2 + s[i]^2)} */ void frb_from_log(frb_signal_t *s); /* Undoes the effect of {frb_to_log}, except for the biasing and abs. Namely, replaces each element {s[i]}, by its exponential {exp(s[i])}. */ void frb_to_log (double eps, frb_signal_t *s) { int nr = s->nel; int k; affirm(eps > 0.0, "bad eps"); for (k = 0; k < nr; k++) { double sk = s->el[k]; sk = 0.5*log(eps*eps + sk*sk); } } void frb_from_log(frb_signal_t *s) { int nr = s->nel; int k; for (k = 0; k < nr; k++) { double *sk = &(s->el[k]); *sk = exp(*sk); } } frb_signal_t **ffta; /* (OUT) Fourier transfs of all {a} shape functions. */ frb_signal_t **pwra; /* (OUT) Unfolded spectra of all {a} shape functions. */ frb_signal_t **fftb; /* (OUT) Fourier transfs of all {b} shape functions. */ frb_signal_t **pwrb; /* (OUT) Unfolded spectra of all {b} shape functions. */ frb_signal_t **fftm; /* (OUT) Fourier transforms of all mean signals {m}. */ frb_signal_t **pwrm; /* (OUT) Unfolded spectra of all mean signals {m}. */ frb_signal_t **fftd; /* (OUT) Fourier transforms of all diff signals {m}. */ frb_signal_t **pwrd; /* (OUT) Unfolded spectra of all diff signals {m}. */ frb_signal_t *tffts; /* Total Fourier transf of {a,b} signals. */ frb_signal_t *tpwrs; /* Total power spectra of {a,b} signals. */ frb_signal_t *tfftm; /* Total Fourier transf of {m} signals. */ frb_signal_t *tpwrm; /* Total power spectra of {m} signals. */ frb_signal_t *tfftd; /* Total Fourier transf of {d} signals. */ frb_signal_t *tpwrd; /* Total power spectra of {d} signals. */ bool_t logScale; /* TRUE accumulates logs of power spectra. */ double eps; /* Noise bias, used when {logScale=TRUE}. */ frb_signal_t **a; /* (OUT) All {a} shape functions. */ frb_signal_t **b; /* (OUT) All {b} shape functions. */ frb_signal_t **m; /* (OUT) All mean signals {m}. */ frb_signal_t **d; /* (OUT) All diff signals {d}. */ AccumSignal(&fft, tfft, FALSE, 0.0); AccumSignal(&pwr, tpwr, logScale, eps); bool_t logScale, /* TRUE accumulates logs of power spectra. */ double eps /* Noise bias, used when {logScale=TRUE}. */ frb_signal_t tffts = frb_signal_new(ns); frb_zero_signal(&tffts); frb_signal_t tpwrs = frb_signal_new(ns); frb_zero_signal(&tpwrs); frb_signal_t tfftm = frb_signal_new(ns); frb_zero_signal(&tfftm); frb_signal_t tpwrm = frb_signal_new(ns); frb_zero_signal(&tpwrm); frb_signal_t tfftd = frb_signal_new(ns); frb_zero_signal(&tfftd); frb_signal_t tpwrd = frb_signal_new(ns); frb_zero_signal(&tpwrd); for (icand = 0; iCand < ncd; iCand++) { frb_signals_from_candidate ( iCand, cvd, ncv, dir, oldCmt, &tffts, &tpwrs, &tfftm, &tpwrm, &tfftd, &tpwrd, logScale, eps ); } ScaleSignal(&tffts, 1.0/(2*ncd)); ScaleSignal(&tpwrs, 1.0/(2*ncd-1)); ScaleSignal(&tfftm, 1.0/(ncd)); ScaleSignal(&tpwrm, 1.0/(ncd-1)); ScaleSignal(&tfftd, 1.0/(ncd)); ScaleSignal(&tpwrd, 1.0/(ncd-1)); int nf = (ns+1)/2; for (k = 0; k < ns; k++) { int freq = (k > frb_dmax ? k - ns : k); double lambda = L/freq; /* Average and variance of the original signals {\bar A_k, \hat A_k}: */ double avgs, vars; frb_compute_avg_var(2*ncd, tffts.el[k], tpwrs.el[k], &avgs, &vars, logScale); /* Average and variance of the mean signals {\bar M, \hat M}: */ double avgm, varm; frb_compute_avg_var(ncd, tfftm.el[k], tpwrm.el[k], &avgm, &varm, logScale); /* Average and variance of the difference signals {\bar D, \hat D}: */ double avgd, vard; frb_compute_avg_var(ncd, tfftd.el[k], tpwrd.el[k], &avgd, &vard, logScale); } free(fft.el); free(pwr.el); argparser_get_keyword(pp, "-tag"); o->tag = argparser_get_next(pp); argparser_get_keyword(pp, "-nFreqs"); o->nFreqs = argparser_get_next_int(pp, 1); argparser_get_keyword(pp, "-nFiles"); o->nFiles = argparser_get_next_int(pp, 1); argparser_get_keyword(pp, "-curvePrefix"); o->curvePrefix = argparser_get_next(pp); if (argparser_keyword_present("-curveDir")) { o->curveDir = argparser_get_next(pp); } else { o->curveDir = "."; } if (argparser_keyword_present("-extension")) { o->extension = argparser_get_next(pp); } else { o->extension = ".flc"; } argparser_get_keyword(pp, "-output"); o->output = argparser_get_next(pp); /* Matching parameters: */ argparser_get_keyword(pp, "-maxShift"); o->maxShift = pp.getNextLongReal(); /* Segment extraction parameters: */ argparser_get_keyword(pp, "-nSteps"); (o->nSamples-1)= argparser_get_next_int(pp, 2,4096); argparser_get_keyword(pp, "-inName"); o->inName = argparser_get_next(pp); argparser_get_keyword(pp, "-outName"); o->outName = argparser_get_next(pp); argparser_get_keyword(pp, "-name"); o->name = argparser_get_next(pp); argparser_get_keyword(pp, "-ixMin"); o->ixMin = argparser_get_next_int(pp, FIRST(INT),LAST(INT)); argparser_get_keyword(pp, "-ixMax"); o->ixMax = argparser_get_next_int(pp, o->iMin,o->iMin+30); argparser_get_keyword(pp, "-iyyMin"); o->iyMin = argparser_get_next_int(pp, FIRST(INT),LAST(INT)); argparser_get_keyword(pp, "-iMax"); o->iyMax = argparser_get_next_int(pp, o->iMin,o->iMin+30); if (argparser_keyword_present("-avgFile")) { o->avgFile = argparser_get_next(pp); } else { o->avgFile = ""; } argparser_skip_parsed(); int nFiles = argc - pp.next; if (nFiles < 2) { argparser_error("needs at least two files"); } o->inFile = malloc(nFiles * sizeof(char *)); for (i = 0; i < nFiles; i++) { char *f = argparser_get_next(pp); o->inFile[i] = f; if ((strlen(f) > 1) && (f[0] == '-')) { argparser_error(txtcat ("bad file name :", f)); } } typedef struct frb_freq_data_t { double freq; double lambda; double rms; } FreqPower = ARRAY OF frb_freq_data_t; typedef struct frb_mismatch_options_t { char * outName; /* Output file name (without extension). */ int nFreqs; /* Number of frequences in output file */; } void frb_compute_band_spectrum(int argc, char**argv) /* Writes the Fourier spectrum (coefs and power) of a uniformly sampled signal. */ /* This program reads a sequence of samples of a real-valued signal f(t). It is assumed that the sampling interval {dt} is constant, and the signal has period {ns*dt} where {ns} is a power of two. The number of samples must be {ns} or {ns+1}; if the latter, the first and last samples must be equal. The program computes the Fourier transform of the sampled signal in the ``smart'' basis {basis(ns,k)[i] = sqrt(2/ns) * cos((2*Pi*k*i / ns) + Pi/4)} and writes it to disk in two formats: basic coefficients {F[k]} (extension {.fft}), and power {P[k]} for each frequency {k} (extension {.fpw}). The power is defined as {P[0] = F[0]^2} {P[k] = F[k]^2 + F[ns-k]^2} for {k} in {1..ns/2-1}, {P[ns/2] = F[ns/2]^2} */ { int k; nFreqs = o->nFreqs; nSpectra = NEW(REF doubleS, nFreqs)^; accPower = NEW(REF doubleS, nFreqs)^; freq = NEW(REF doubleS, nFreqs)^; lambda = NEW(REF doubleS, nFreqs)^; accMean = NEW(REF doubleS,nFreqs)^ for (k = 0; k <= nFreqs-1; k++) { nSpectra[k] = 0.0; accPower[k] = 0.0; freq[k] = 0.0; lambda[k] = 0.0; } LOOP Lex.Skip(stdin); if (Rd.EOF(stdin)) { EXIT; } frb_read_and_accumulate_freq_spectrum (stdin, nSpectra, accPower, freq, lambda, nFreqs); } k= 0; while ( (k < nFreqs)) { if (nSpectra[k] > 0.0) { accMean[k] = accPower[k]/nSpectra[k]; k = k+1; } } frb_write_mean_freq_power(o->outName, accMean, freq, lambda); } void frb_read_and_accumulate_freq_spectrum( Rd.T rd; doubleS nSpectra; doubleS accPower; doubleS freq; doubleS lambda; int nFreqs) = { { ?? s = frb_read_freq_power(rd, nFreqs)^; n = s->nel ; for (k = 0; k <= n-1; k++) { nSpectra[k] = nSpectra[k] + 1.0; accPower[k]= accPower[k] + s[k].rms*s[k].rms; lambda[k] = s[k].lambda; freq[k] = s[k].freq; } } } REF FreqPower frb_read_freq_power (Rd.T rd; int nFreqs) = int k; { { ?? /* rd = open_read(fileName), */ su = NEW(REF FreqPower, nFreqs); s = su^ ; EVAL FileFmt.ReadComment(rd, '#'); /* fprintf(stderr, "reading file ...\n"); */ /* fprintf(stderr, cmt); fprintf(stderr, "\n"); */ k= 0; FGet.Skip(rd); for (k = 0; k <= nFreqs-1; k++) { s[k].freq = FGet.LongReal(rd); /* fprintf(stderr, Fmt.LongReal(s[k].freq) & " " & "\n"); */ s[k].lambda = FGet.LongReal(rd); /* fprintf(stderr, Fmt.LongReal(s[k].lambda) & " "); */ s[k].rms = FGet.LongReal(rd); /* fprintf(stderr, Fmt.LongReal(s[k].rms) & "\n"); */ FGet.Skip(rd); } return su; } } void frb_write_mean_freq_power( char * outName; doubleS *accMean; doubleS *freq; doubleS *lambda; )= { { ?? outCandName = outName & ".fpw"; wr = open_write(outCandName) ; fprintf(stderr, "writing file " & outCandName & " ...\n"); fprintf(wr, "# freq lambda mean \n"); for (k = 0; k <= LAST(accMean); k++) { if (lambda[k] >= 0.0) { FPut.LongReal(wr, freq[k], 2, Fmt.Style.Fix); fprintf(wr, " "); FPut.LongReal(wr, lambda[k], 2, Fmt.Style.Fix); fprintf(wr, " "); FPut.LongReal(wr, sqrt(accMean[k]), 4, Fmt.Style.Fix); fprintf(wr, "\n"); } } Wr.Flush(wr); } } #define MaxBands = 256; typedef BandData = RECORD double lo; double hi; double rms; } Spectrum = ARRAY OF BandData; frb_mismatch_options_t = PZMismatch.frb_options_t;struct frb_options_t { char * inName; /* Prefix for input file names. */ int nFiles; /* Number of files to read. */ char * outName; /* Prefix for output file names. */ char * tag; /* Suffix for input and output file names. */; } void frb_compute_band_spectrum_2(int argc, char**argv) /* This program reads one or more band-power spectrum files with names of the form `xxx-iiiiii-yyy.bpw' where `xxx' is {o->inName}, `yyy' is {o->tag}, and `iiiiii' is a six-digit serial number. Each file must contain sero or more records of the form "LAMBDALO LAMBDAHI RMS" where {RMS} is the square root of the total power in the wavelength range {[LAMBDALO _ LAMBDAHI)}. All spectra must contain the same bands in the same order, except that some bands at the end may be missing in some files. The program computes the power in each band, averaged over all spectra that include that band, and writes the square root of the result as a single file `zzz-mean-yyy.bpw' where `zzz' is {o->outName}. */ { int k; nFiles = o->nFiles; tag = o->tag; nSpectra = NEW(REF doubleS, MaxBands)^; accPower = NEW(REF doubleS, MaxBands)^; lambdaLo = NEW(REF doubleS, MaxBands)^; lambdaHi = NEW(REF doubleS, MaxBands)^; accMean = NEW(REF doubleS,MaxBands)^ for (k = 0; k <= MaxBands-1; k++) { nSpectra[k] = 0.0; accPower[k] = 0.0; lambdaLo[k] = 0.0; lambdaHi[k] = 0.0; } for (i = 0; i <= nFiles-1; i++) { inFile = o->inName & "-" & Fmt.Pad(Fmt.Int(i), 6, '0') & "-" & tag & ".bpw" frb_read_and_accumulate_spectrum(inFile, nSpectra, accPower, lambdaLo, lambdaHi); } k= 0; while ( nSpectra[k] > 0.0) { accMean[k] = accPower[k]/nSpectra[k]; k = k+1; } outCandName = o->outName & "-mean-" & tag & ".bpw" frb_write_mean_band_power(outCandName, accMean, lambdaLo, lambdaHi); } void frb_read_and_accumulate_spectrum( char * fileName; doubleS nSpectra; doubleS accPower; doubleS lambdaLo; doubleS lambdaHi; ) = { { ?? s = frb_read_band_power(fileName)^; n = s->nel ; for (k = 0; k <= n-1; k++) { nSpectra[k] = nSpectra[k] + 1.0; accPower[k]= accPower[k] + s[k].rms*s[k].rms; if (lambdaHi[k] == 0.0) { lambdaLo[k] = s[k].lo; lambdaHi[k] = s[k].hi; affirm(lambdaLo[k] < lambdaHi[k] , "bug"); affirm(lambdaLo[k] >= 0.0 , "bug");; } else { affirm(lambdaLo[k] == s[k].lo , "bug"); affirm(lambdaHi[k] == s[k].hi , "bug");; } } } } REF Spectrum frb_read_band_power (char * fileName) = int k; { { ?? rd = open_read(fileName); s = NEW(REF Spectrum, MaxBands)^; cmt = FileFmt.ReadComment(rd, '#') ; fprintf(stderr, "reading file " & fileName & " ...\n"); fprintf(stderr, cmt); fprintf(stderr, "\n"); k= 0; REPEAT FGet.Skip(rd); s[k].lo = FGet.LongReal(rd); fprintf(stderr, "Low lambda: " & Fmt.LongReal(s[k].lo) & "\n"); s[k].hi = FGet.LongReal(rd); fprintf(stderr, "Hi lambda: " & Fmt.LongReal(s[k].hi) & "\n"); s[k].rms = FGet.LongReal(rd); fprintf(stderr, "sqrt(power) :" & Fmt.LongReal(s[k].rms) & "\n"); FGet.Skip(rd); k = k + 1 UNTIL Rd.EOF(rd); { ?? w = NEW(REF Spectrum, k) ; w^ = SUBARRAY(s, 0, k); return w; } } } void frb_write_mean_band_power( char * outCandName; doubleS *accMean; doubleS *lambdaLo; doubleS *lambdaHi; )= { { ?? wr = open_write(outCandName) ; fprintf(stderr, "writing file " & outCandName & " ...\n"); fprintf(wr, "# lambdaLo lambdaHi mean \n"); for (k = 0; k <= LAST(accMean); k++) { if (lambdaHi[k] > 0.0) { FPut.LongReal(wr, lambdaLo[k], 2, Fmt.Style.Fix); fprintf(wr, " "); FPut.LongReal(wr, lambdaHi[k], 2, Fmt.Style.Fix); fprintf(wr, " "); FPut.LongReal(wr, sqrt(accMean[k]), 4, Fmt.Style.Fix); fprintf(wr, "\n"); } } Wr.Flush(wr); } } void frb_compute_packed_stats(int argc, char**argv) { double v; nFiles = o->inFile^->nel; frb_signal_t rchs[o->nFiles]; frb_signal_t avg = frb_avg_signals(rchs, o->nFiles, o->nSamples); frb_signal_t var = frb_var_signals(rchs, &avg, o->nFiles, o->nSamples); frb_curve_t res = frb_pack_stats(&avg, &var); frb_signal_t *rchs = frb_read_n_signals(o->nFiles, o->nSamples, o->inFile, log); if (o->log) { frb_from_log(&res); } frb_curve_write ( open_write(o->outCandName), cmt = "average and confidence interval of the following files:\n" & frb_make_file_name_list(o->inFile^), c = res, unit = frb_select_unit(res, FRB_MAX_SCALED) ); } char *frb_make_file_name_list(char **file, int nFiles) /* Makes a single string from several file names, one per line. */ { int i; /* Compute total length of all strings: */ int totlen = 0; for (i = 0; i < nFiles; i++) { totlen += strlen(file[i]); } /* Concatenate them: */ char *res = malloc(totlen + 3*nFiles + 1); char *p = res; for (i = 0; i < nFiles; i++) { int size = sprintf(p, " %s\n", file[i]); p += size; } return res; } (frb_match_packed_t){0, (i2_t){{0,0}}, (i2_t){{0,0}}, NULL} frb_signal_t frb_signal_trim ( frb_signal_t *c, int start, int length ) { frb_signal_t x = double_vec_new(length); frb_signal_do_trim(c, start, length, &x); return x; } void frb_signal_do_trim ( frb_signal_t *c, int start, int length, frb_signal_t *x ) { int n = c->nel; int ini = imod(start, n); int i, j; for (i= 0, j = ini; i < length; i++, j++) { if (j > n) { j -= n; } x->el[i] = c->el[j]; } } double frb_signal_adjust_unit ( double givenUnit, frb_signal_t *c ) { int n = c->nel; double avg = 0, dev = 0, big = 0; /* Compute average and max-abs value: */ if (n > 0 ) { int i; for (i = 0; i < n; i++) { double ci = c->el[i]; avg = avg + ci; if (abs(ci) > big) { big = abs(ci); } } avg = avg/(double)n; } /* Compute standard deviation: */ if (n > 1) { int i; for (i = 0; i < n; i++) { double d = c->el[i] - avg; dev += d*d; } dev = sqrt(dev/((double)n-1)); } return frb_adjust_unit(givenUnit, dev, big); } int band, char *extension void frb_curve_fourier_transform ( frb_curve_t *c, frb_curve_t *f ) { int nc = (c->nel); affirm(f->nel == nc, "wrong size"); frb_signal_t cr = double_vec_new(nc); frb_signal_t fr = double_vec_new(nc); int i, k; for (k = 0; k < 3; k++) { for (i = 0; i < nc; i++) { cr.el[i] = c->el[i].c[k]; } frb_fourier_transform(&cr, &fr); for (i = 0; i < nc; i++) { f->el[i].c[k] = fr.el[i]; } } } void frb_curve_inv_fourier_transform ( frb_curve_t *f, frb_curve_t *c ) { int nc = (c->nel); affirm(f->nel == nc, "wrong size"); frb_signal_t cr = double_vec_new(nc); frb_signal_t fr = double_vec_new(nc); int i, k; for (k = 0; k < 3; k++) { for (i = 0; i < nc; i++) { fr.el[i] = f->el[i].c[k]; } frb_fourier_transform(&fr, &cr); for (i = 0; i < nc; i++) { c->el[i].c[k] = cr.el[i]; } } } void frb_compute_avg_var ( int n, double sum, double sum2, double *avg, double *var, bool_t logScale ) { *avg = sum/n; *var = sum2/n; if (logScale) { *var = exp(*var); } } void frb_compute_avg_var(int n, double sum, double sum2, double *avg, double *var, bool_t logScale); /* Computes the average and variance of a list of values, given the number {n} of values, their sum {sum} and the sum of their squares {sum2}. The variance is computed assuming that the mean value is 0.0. If {logScale=TRUE}, assumes that {sum2} is the sum of the logs of the samples squared, and returns the {exp} of the result. */