home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: Science / Science.zip / gmt_os2.zip / src / spectrum1d.c < prev    next >
C/C++ Source or Header  |  1995-03-13  |  17KB  |  556 lines
  1. /*--------------------------------------------------------------------
  2.  *    The GMT-system: @(#)spectrum1d.c    2.12  3/13/95    
  3.  *
  4.  *    Copyright (c) 1991-1995 by P. Wessel and W. H. F. Smith
  5.  *    See README file for copying and redistribution conditions.
  6.  *--------------------------------------------------------------------*/
  7. /*
  8.  *
  9.  *    spectrum1d <file> [-C] [-D<dt>] [-N<name_stem>] -S<segment_size> [-V] [-W]
  10.  *
  11.  *    Compute auto and cross spectra using Welch's method of
  12.  *    multiple overlapped windows.  Find 1 standard error bars
  13.  *    following expressions in Bendat & Piersol.
  14.  *
  15.  *    -D<dt>    set delta_time, the sampling of the timeseries.
  16.  *        [Default <dt> = 1.0]
  17.  *
  18.  *    -N<name_stem>    name stem for filenames.  Files will be
  19.  *        created called <name_stem>.xpower, etc.
  20.  *        [Default <name_stem> = "spectrum"]
  21.  *
  22.  *    -S<segment_size>    give an integer radix 2 window width.
  23.  *        Total timeseries will be split into pieces of
  24.  *        length <segment_size>.  Std errors in spectra are
  25.  *        approximately 1/sqrt(n_data/segment_size).
  26.  *
  27.  *    -V    Verbose operation; write info to stderr.
  28.  *        [Default is silent]
  29.  *
  30.  *    -W    write Wavelength in col 1 instead of frequency.
  31.  *        [Default writes frequency in cycles/dt_units]
  32.  *
  33.  *    -C    input has 2 cols, X(t) Y(t); do cross-spectra.
  34.  *        [Default is one column; do X power spectrum only]
  35.  *
  36.  *    Author:        W. H. F. Smith
  37.  *    Date:        11 April 1991-1995
  38.  *    Revised:    5 June 1991-1995 to add W and N options in prep for
  39.  *            GMT v2.0 release.
  40.  *    References:    Julius S. Bendat & Allan G. Piersol,
  41.  *            "Random Data", 2nd revised edition, 566pp.,
  42.  *            1986, John Wiley & Sons, New York.
  43.  *
  44.  *            Peter D. Welch, "The use of Fast Fourier
  45.  *            Transform for the estimation of power spectra:
  46.  *            a method based on time averaging over short,
  47.  *            modified periodograms", IEEE Transactions on
  48.  *            Audio and Electroacoustics, Vol AU-15, No 2,
  49.  *            June, 1967.
  50.  */
  51.  
  52. #include "gmt.h"
  53.  
  54. float    *datac, *x, *y;
  55.  
  56. struct    SPEC {
  57.         double    xpow;    /* PSD in X(t)  */
  58.         double    ypow;    /* PSD in Y(t)  */
  59.         double    gain;    /* Amplitude increase X->Y in Optimal Response Function  */
  60.         double    phase;    /* Phase of Y w.r.t. X in Optimal Response Function  */
  61.         double    coh;    /* (squared) Coherence of Y and X; SNR of Y = coh/(1-coh)  */
  62.         double    radmit;    /* Real part of Admittance; used e.g., for gravity/topography  */
  63. }    *spec;
  64.  
  65. double    dt = 1.0, x_variance, y_variance, d_n_windows, y_pow;
  66.  
  67. int    y_given = FALSE, write_wavelength = FALSE;
  68. int    n_spec,  window = 0, window_2, n_data;
  69.  
  70. char    namestem[80];
  71.  
  72. FILE    *fp = NULL;
  73.  
  74. void    alloc_arrays();
  75. void    read_data();
  76. void    compute_spectra();
  77. void    detrend_and_hanning();
  78. void    write_output();
  79. void    free_space();
  80.  
  81. main (argc, argv)
  82. int    argc;
  83. char    **argv;
  84. {
  85.  
  86.     int    i, window_test, error = FALSE;
  87.     
  88.     argc = gmt_begin (argc, argv);
  89.     
  90.     sprintf(namestem, "spectrum\0");
  91.  
  92.     for (i = 1; i < argc; i++) {
  93.         if (argv[i][0] == '-') {
  94.             switch (argv[i][1]) {
  95.         
  96.                 /* Common parameters */
  97.             
  98.                 case 'H':
  99.                 case 'V':
  100.                 case '\0':
  101.                     error += get_common_args (argv[i], 0, 0, 0, 0);
  102.                     break;
  103.                 
  104.                 /* Supplemental parameters */
  105.                 
  106.                 case 'C':
  107.                     y_given = TRUE;
  108.                     break;
  109.                 case 'D':
  110.                     dt = atof(&argv[i][2]);
  111.                     break;
  112.                 case 'N':
  113.                     strcpy(namestem, &argv[i][2]);
  114.                     if (!namestem[0]) error = TRUE;
  115.                     break;
  116.                 case 'S':
  117.                     window = atoi(&argv[i][2]);
  118.                     window_test = 2;
  119.                     while (window_test < window) {
  120.                         window_test += window_test;
  121.                     }
  122.                     break;
  123.                 case 'W':
  124.                     write_wavelength = TRUE;
  125.                     break;
  126.                 default:
  127.                     error = TRUE;
  128.                     gmt_default_error (argv[i][1]);
  129.                     break;
  130.             }
  131.         }
  132.         else {
  133.             if ((fp = fopen(argv[i], "r")) == NULL) {
  134.                 fprintf(stderr,"spectrum1d:  Cannot open r %s\n", argv[i]);
  135.                 error = TRUE;
  136.             }
  137.         }
  138.     }
  139.  
  140.     if (argc == 1 || gmt_quick) {    /* Display usage */
  141.         fprintf (stderr, "spectrum1d %s - compute auto- [and cross- ] spectra from one [or two] timeseries\n\n", GMT_VERSION);
  142.         fprintf(stderr,"usage:  spectrum1d -S<segment_size> [-C] [-D<dt>] [-H] [-N<name_stem>] [-V] [-W]\n");
  143.         
  144.         if (gmt_quick) exit (-1);
  145.         
  146.         fprintf(stderr,"\t-S<segment_size>  Use data subsets of <segment_size> elements.\n");
  147.         fprintf(stderr,"\t\t<segment_size> must be radix 2; std. err. = 1/sqrt(n_data/segment_size).\n");
  148.         fprintf(stderr,"\tOptions:\n");
  149.         fprintf(stderr,"\t-C  2 column X(t),Y(t) input; estimate Cross-spectra [Default 1 col, X power only].\n\n");
  150.         fprintf(stderr,"\t-D<dt>  set delta_time sampling interval of data [Default = 1.0].\n");
  151.         explain_option ('H');
  152.         fprintf(stderr,"\t-N<name_stem>  supply name stem for files [Default = 'spectrum'].\n");
  153.         fprintf(stderr,"\t\tOutput files will be named <name_stem>.xpower, etc.\n");
  154.         explain_option ('V');
  155.         fprintf(stderr,"\t-W  write Wavelength of spectral estimate in col 1 [Default = frequency].\n");
  156.         exit(-1);
  157.     }
  158.  
  159.     if (window == 0) {
  160.         fprintf (stderr, "%s: GMT SYNTAX ERROR -S option: segment size must be positive\n", gmt_program);
  161.         error++;
  162.     }
  163.     if (window_test != window) {
  164.         fprintf (stderr, "%s: GMT SYNTAX ERROR -S option: Segment size not radix 2.  Try %d or %d\n", gmt_program,
  165.             (window_test/2), window_test);
  166.         error++;
  167.     }
  168.     if (dt <= 0.0) {
  169.         fprintf (stderr, "%s: GMT SYNTAX ERROR -D option: Sampling interval must be positive\n", gmt_program);
  170.         error++;
  171.     }
  172.  
  173.     if (error) exit (-1);
  174.     
  175.     alloc_arrays();
  176.  
  177.     read_data();
  178.  
  179.     compute_spectra();
  180.  
  181.     write_output();
  182.  
  183.     free_space();
  184.  
  185.     gmt_end (argc, argv);
  186. }
  187.  
  188. void    alloc_arrays()
  189. {
  190.     n_spec = window/2;    /* This means we skip zero frequency; data are detrended  */
  191.     window_2 = 2 * window;        /* This is for complex array stuff  */
  192.  
  193.     spec = (struct SPEC *) memory (CNULL, n_spec, sizeof(struct SPEC), "spectrum1d");
  194.     datac = (float *) memory (CNULL, window_2, sizeof(float), "spectrum1d");
  195. }
  196.  
  197. void    read_data()
  198. {
  199.     int    i, n_alloc;
  200.     char    buffer[512];
  201.     double    xx, yy;
  202.  
  203.     if (fp == NULL) fp = stdin;
  204.     n_alloc = GMT_CHUNK;
  205.     n_data = 0;
  206.     x = (float *) memory (CNULL, n_alloc, sizeof(float), "spectrum1d");
  207.     if (y_given) y = (float *)memory (CNULL, n_alloc, sizeof(float), "spectrum1d");
  208.  
  209.     if (gmtdefs.io_header) for (i = 0; i < gmtdefs.n_header_recs; i++) fgets (buffer, 512, fp);
  210.  
  211.         while ( (fgets (buffer, 512, fp) ) != NULL) {
  212.  
  213.                 if ((sscanf(buffer, "%lf %lf", &xx, &yy)) != (1 + y_given)) {
  214.                         fprintf(stderr,"spectrum1d:  Cannot read line %d\n", n_data+1);
  215.                         continue;
  216.                 }
  217.                 x[n_data] = xx;
  218.         if (y_given) y[n_data] = yy;
  219.                 n_data++;
  220.  
  221.                 if (n_data == n_alloc) {
  222.                         n_alloc += GMT_CHUNK;
  223.             x = (float *) memory ((char *)x, n_alloc, sizeof(float), "spectrum1d");
  224.             if (y_given) y = (float *)memory ((char *)y, n_alloc, sizeof(float), "spectrum1d");
  225.                 }
  226.         }
  227.         if (fp != stdin) fclose(fp);
  228.     x = (float *) memory ((char *)x, n_data, sizeof(float), "spectrum1d");
  229.     if (y_given) y = (float *)memory ((char *)y, n_data, sizeof(float), "spectrum1d");
  230.     if (gmtdefs.verbose) fprintf(stderr,"Read %d data points.\n", n_data);
  231. }
  232.  
  233. void    compute_spectra()
  234. {
  235.     int    n_windows, w, i, t_start, t_stop, t, f;
  236.     int    narray, ndim = 1, ksign = -1, iform = 1;
  237.     float    work = 0.0;
  238.     double    dw, spec_scale, x_varp, y_varp, one_on_nw, co_quad;
  239.     double    xreal, ximag, yreal, yimag, xpower, ypower, co_spec, quad_spec;
  240.  
  241.     narray = window;
  242.  
  243.     /* Scale factor for spectral estimates should be 1/4 of amount given in
  244.         Bendat & Piersol eqn 11-102 because I compute 2 * fft in my
  245.         one-sided code below.  However, tests show that I need 1/8 of
  246.         their equation to match variances approximately:  */
  247.  
  248.     /* This used to read:  spec_scale = 0.5 / (window_2 * dt);  */
  249.     spec_scale = dt / (window_2);
  250.  
  251.     d_n_windows = (double)n_data / (double)window;
  252.  
  253.     n_windows = rint(2.0 * d_n_windows) - 1;
  254.     one_on_nw = 1.0 / n_windows;
  255.     dw = (double)(n_data - window) / (double)(n_windows - 1);
  256.  
  257.     for (w = 0; w < n_windows; w++) {
  258.         t_start = floor (0.5 + w * dw);
  259.         t_stop = t_start + window;
  260.         if (y_given) {
  261.             for (t = t_start, i = 0; t < t_stop; t++, i+=2) {
  262.                 datac[i] = x[t];
  263.                 datac[i+1] = y[t];
  264.             }
  265.         }
  266.         else {
  267.             for (t = t_start, i = 0; t < t_stop; t++, i+=2) {
  268.                 datac[i] = x[t];
  269.                 datac[i+1] = 0.0;
  270.             }
  271.         }
  272.  
  273.         detrend_and_hanning();
  274.  
  275.         fourt_(datac, &narray, &ndim, &ksign, &iform, &work);
  276.  
  277.         /* Get one-sided estimates:  */
  278.  
  279.         x_varp = spec_scale * (datac[0] * datac[0]);
  280.         if (y_given) {
  281.             y_varp = spec_scale * (datac[1] * datac[1]);
  282.             for (i = 0, f = 2; i < n_spec; i++, f+=2) {
  283.                 xreal = (i == n_spec - 1) ? datac[f] : datac[f] + datac[window_2 - f];
  284.                 ximag = (i == n_spec - 1) ? 0.0 : datac[f+1] - datac[window_2 - f + 1];
  285.                 yreal = (i == n_spec - 1) ? datac[f+1] : datac[f+1] + datac[window_2 - f + 1];
  286.                 yimag = (i == n_spec - 1) ? 0.0 : datac[window_2 - f] - datac[f];
  287.                 xpower = spec_scale * (xreal * xreal + ximag * ximag);
  288.                 ypower = spec_scale * (yreal * yreal + yimag * yimag);
  289.                 co_spec = spec_scale * (xreal * yreal + ximag * yimag);
  290.                 quad_spec = spec_scale * (ximag * yreal - yimag * xreal);
  291.  
  292.                 x_varp += xpower;
  293.                 y_varp += ypower;
  294.                 spec[i].xpow += xpower;
  295.                 spec[i].ypow += ypower;
  296.                 /* Temporarily store co-spec in gain:  */
  297.                 spec[i].gain += co_spec;
  298.                 /* Temporarily store quad-spec in phase:  */
  299.                 spec[i].phase += quad_spec;
  300.             }
  301.             x_varp *= (dt/n_spec);
  302.             y_varp *= (dt/n_spec);
  303.         }
  304.         else {
  305.             for (i = 0, f = 2; i < n_spec; i++, f+=2) {
  306.                 xreal = datac[f] + datac[window_2 - f];
  307.                 ximag = datac[f+1] - datac[window_2 - f + 1];
  308.                 xpower = spec_scale * (xreal * xreal + ximag * ximag);
  309.                 x_varp += xpower;
  310.                 spec[i].xpow += xpower;
  311.             }
  312.             x_varp *= (dt/n_spec);
  313.         }
  314.  
  315.         if (gmtdefs.verbose) {
  316.             y_pow = (y_given) ? y_variance/y_varp : 0.0;
  317.             fprintf(stderr,"Window %d from %d to %d\n", w, t_start, t_stop);
  318.             fprintf(stderr,"X var:  %.8lg  X pow:  %.8lg  ratio:  %.8lg  Y var:  %.8lg  Y pow:  %.8lg  ratio:  %.8lg\n",
  319.                 x_variance, x_varp, (x_variance/x_varp), y_variance, y_varp, y_pow);
  320.         }
  321.     }
  322.     /* Now we can divide by n_windows for the ensemble average.
  323.         The cross spectral stuff needs to be computed:  */
  324.  
  325.     if (y_given ) {
  326.         for (i = 0; i < n_spec; i++) {
  327.             spec[i].xpow *= one_on_nw;
  328.             spec[i].ypow *= one_on_nw;
  329.             co_spec = spec[i].gain * one_on_nw;
  330.             quad_spec = spec[i].phase * one_on_nw;
  331.             spec[i].phase = d_atan2(quad_spec, co_spec);
  332.             co_quad = co_spec * co_spec + quad_spec * quad_spec;
  333.             spec[i].coh = co_quad / (spec[i].xpow * spec[i].ypow);
  334.             spec[i].gain = sqrt(co_quad) / spec[i].xpow;
  335.             spec[i].radmit = co_spec / spec[i].xpow;
  336.         }
  337.     }
  338.     else {
  339.         for (i = 0; i < n_spec; i++) {
  340.             spec[i].xpow *= one_on_nw;
  341.         }
  342.     }
  343. }
  344.  
  345. void    write_output()
  346. {
  347.     int    i;
  348.     double    f, delta_f, eps_pow, eps, cop, nop;
  349.     char    fname[120], format[50];
  350.     FILE    *fpout;
  351.  
  352.     sprintf (format, "%s\t%s\t%s\n\0", gmtdefs.d_format, gmtdefs.d_format, gmtdefs.d_format);
  353.     
  354.     delta_f = 1.0 / (window * dt);
  355.     eps_pow = 1.0 / sqrt(d_n_windows);    /* Multiplicative error bars for power spectra  */
  356.  
  357.     /* First write x power  */
  358.  
  359.     sprintf(fname, "%s.xpower\0", namestem);
  360.     if ( (fpout = fopen(fname, "w")) == NULL) {
  361.         fprintf(stderr,"spectrum1d:  Cannot open w %s\n", fname);
  362.         exit(-1);
  363.     }
  364.     if (gmtdefs.verbose) fprintf(stderr,"spectrum1d:  Writing %s\n", fname);
  365.     for (i = 0; i < n_spec; i++) {
  366.         f = (i + 1) * delta_f;
  367.         if (write_wavelength) f = 1.0 / f;
  368.         fprintf(fpout, format, f, spec[i].xpow, eps_pow * spec[i].xpow);
  369.     }
  370.     fclose(fpout);
  371.  
  372.     if (y_given) {
  373.         /* Write y power  */
  374.         sprintf(fname, "%s.ypower\0", namestem);
  375.         if ( (fpout = fopen(fname, "w")) == NULL) {
  376.             fprintf(stderr,"spectrum1d:  Cannot open w %s\n", fname);
  377.             exit(-1);
  378.         }
  379.         if (gmtdefs.verbose) fprintf(stderr,"spectrum1d:  Writing %s\n", fname);
  380.         for (i = 0; i < n_spec; i++) {
  381.             f = (i + 1) * delta_f;
  382.             if (write_wavelength) f = 1.0 / f;
  383.             fprintf(fpout, format, f, spec[i].ypow, eps_pow * spec[i].ypow);
  384.         }
  385.         fclose(fpout);
  386.  
  387.         /* Write Coherent Output power  */
  388.         sprintf(fname, "%s.cpower\0", namestem);
  389.         if ( (fpout = fopen(fname, "w")) == NULL) {
  390.             fprintf(stderr,"spectrum1d:  Cannot open w %s\n", fname);
  391.             exit(-1);
  392.         }
  393.         if (gmtdefs.verbose) fprintf(stderr,"spectrum1d:  Writing %s\n", fname);
  394.         for (i = 0; i < n_spec; i++) {
  395.             f = (i + 1) * delta_f;
  396.             if (write_wavelength) f = 1.0 / f;
  397.             cop = spec[i].ypow * spec[i].coh;
  398.             eps = eps_pow * sqrt( (2.0 - spec[i].coh) / spec[i].coh);
  399.             fprintf(fpout, format, f, cop, eps*cop);
  400.         }
  401.         fclose(fpout);
  402.  
  403.         /* Write Noise Output power  */
  404.         sprintf(fname, "%s.npower\0", namestem);
  405.         if ( (fpout = fopen(fname, "w")) == NULL) {
  406.             fprintf(stderr,"spectrum1d:  Cannot open w %s\n", fname);
  407.             exit(-1);
  408.         }
  409.         if (gmtdefs.verbose) fprintf(stderr,"spectrum1d:  Writing %s\n", fname);
  410.         for (i = 0; i < n_spec; i++) {
  411.             f = (i + 1) * delta_f;
  412.             if (write_wavelength) f = 1.0 / f;
  413.             nop = spec[i].ypow * (1.0 - spec[i].coh);
  414.             fprintf(fpout, format, f, nop, eps_pow*nop);
  415.         }
  416.         fclose(fpout);
  417.  
  418.         /* Write Gain spectrum  */
  419.         sprintf(fname, "%s.gain\0", namestem);
  420.         if ( (fpout = fopen(fname, "w")) == NULL) {
  421.             fprintf(stderr,"spectrum1d:  Cannot open w %s\n", fname);
  422.             exit(-1);
  423.         }
  424.         if (gmtdefs.verbose) fprintf(stderr,"spectrum1d:  Writing %s\n", fname);
  425.         for (i = 0; i < n_spec; i++) {
  426.             f = (i + 1) * delta_f;
  427.             if (write_wavelength) f = 1.0 / f;
  428.             eps = eps_pow * sqrt( (1.0 - spec[i].coh) / 2.0 * spec[i].coh);
  429.             fprintf(fpout, format, f, spec[i].gain, eps*spec[i].gain);
  430.         }
  431.         fclose(fpout);
  432.  
  433.         /* Write Real Admittance spectrum  */
  434.         sprintf(fname, "%s.admit\0", namestem);
  435.         if ( (fpout = fopen(fname, "w")) == NULL) {
  436.             fprintf(stderr,"spectrum1d:  Cannot open w %s\n", fname);
  437.             exit(-1);
  438.         }
  439.         if (gmtdefs.verbose) fprintf(stderr,"spectrum1d:  Writing %s\n", fname);
  440.         for (i = 0; i < n_spec; i++) {
  441.             f = (i + 1) * delta_f;
  442.             if (write_wavelength) f = 1.0 / f;
  443.             eps = eps_pow * sqrt( (1.0 - spec[i].coh) / 2.0 * spec[i].coh);
  444.             fprintf(fpout, format, f, spec[i].radmit, fabs (eps*spec[i].radmit));
  445.         }
  446.         fclose(fpout);
  447.  
  448.         /* Write Phase spectrum  */
  449.         sprintf(fname, "%s.phase\0", namestem);
  450.         if ( (fpout = fopen(fname, "w")) == NULL) {
  451.             fprintf(stderr,"spectrum1d:  Cannot open w %s\n", fname);
  452.             exit(-1);
  453.         }
  454.         if (gmtdefs.verbose) fprintf(stderr,"spectrum1d:  Writing %s\n", fname);
  455.         for (i = 0; i < n_spec; i++) {
  456.             f = (i + 1) * delta_f;
  457.             if (write_wavelength) f = 1.0 / f;
  458.             eps = eps_pow * sqrt( (1.0 - spec[i].coh) / 2.0 * spec[i].coh);
  459.             fprintf(fpout, format, f, spec[i].phase, eps);
  460.         }
  461.         fclose(fpout);
  462.  
  463.         /* Write Coherency spectrum  */
  464.         sprintf(fname, "%s.coh\0", namestem);
  465.         if ( (fpout = fopen(fname, "w")) == NULL) {
  466.             fprintf(stderr,"spectrum1d:  Cannot open w %s\n", fname);
  467.             exit(-1);
  468.         }
  469.         if (gmtdefs.verbose) fprintf(stderr,"spectrum1d:  Writing %s\n", fname);
  470.         for (i = 0; i < n_spec; i++) {
  471.             f = (i + 1) * delta_f;
  472.             if (write_wavelength) f = 1.0 / f;
  473.             eps = eps_pow * (1.0 - spec[i].coh) * sqrt(2.0 / spec[i].coh);
  474.             fprintf(fpout, format, f, spec[i].coh, eps * spec[i].coh);
  475.         }
  476.         fclose(fpout);
  477.     }
  478. }
  479.  
  480. void    free_space ()
  481. {
  482.     free ((char *)spec);
  483.     free ((char *)datac);
  484.     free ((char *)x);
  485.     if (y_given) free ((char *)y);
  486. }
  487.  
  488. void    detrend_and_hanning()
  489. {
  490.     int    i, t;
  491.     double    sumx, sumtx, sumy, sumty, sumt2, x_slope, x_mean, y_slope, y_mean;
  492.     double    t_factor, h_period, h_scale, hc, hw, tt;
  493.     sumx = 0.0;
  494.     sumtx = 0.0;
  495.     sumy = 0.0;
  496.     sumty = 0.0;
  497.     sumt2 = 0.0;
  498.     x_variance = 0.0;
  499.     y_variance = 0.0;
  500.     t_factor = 2.0 / (window - 1);
  501.     h_period = M_PI / (double)window;    /* For Hanning window  */
  502.     h_scale = sqrt(8.0/3.0);        /* For Hanning window  */
  503.  
  504.     if (y_given) {
  505.         for (i = 0, t = 0; i < window_2; i+=2, t++) {
  506.             tt = t * t_factor - 1.0;
  507.             sumt2 += (tt * tt);
  508.             sumx += datac[i];
  509.             sumtx += (tt * datac[i]);
  510.             sumy += datac[i+1];
  511.             sumty += (tt * datac[i+1]);
  512.         }
  513.     }
  514.     else {
  515.         for (i = 0, t = 0; i < window_2; i+=2, t++) {
  516.             tt = t * t_factor - 1.0;
  517.             sumt2 += (tt * tt);
  518.             sumx += datac[i];
  519.             sumtx += (tt * datac[i]);
  520.         }
  521.     }
  522.     x_slope = sumtx / sumt2;
  523.     x_mean = sumx / window;
  524.     if (y_given) {
  525.         y_slope = sumty / sumt2;
  526.         y_mean = sumy / window;
  527.         for (i = 0, t = 0; i < window_2; i+=2, t++) {
  528.             hc = cos(t * h_period);
  529.             hw = h_scale * (1.0 - hc * hc);
  530.             tt = t * t_factor - 1.0;
  531.             datac[i] -= (x_mean + tt * x_slope);
  532.             datac[i] *= hw;
  533.             x_variance += (datac[i] * datac[i]);
  534.             datac[i+1] -= (y_mean + tt * y_slope);
  535.             datac[i+1] *= hw;
  536.             y_variance += (datac[i+1] * datac[i+1]);
  537.         }
  538.         x_variance /= window;
  539.         y_variance /= window;
  540.     }
  541.     else {
  542.         for (i = 0, t = 0; i < window_2; i+=2, t++) {
  543.             hc = cos(t * h_period);
  544.             hw = h_scale * (1.0 - hc * hc);
  545.             tt = t * t_factor - 1.0;
  546.             datac[i] -= (x_mean + tt * x_slope);
  547.             datac[i] *= hw;
  548.             x_variance += (datac[i] * datac[i]);
  549.         }
  550.         x_variance /= window;
  551.     }
  552. }
  553.  
  554.  
  555.  
  556.