var/external/wavelet/awaprogs/chap05/cd.c

   1 /*
   2  * Convolution-decimation functions and their adjoints.
   3  *
   4  * Copyright (C) 1991--94 Wickerhauser Consulting.  All Rights Reserved.
   5  * May be freely copied for noncommercial use.  See
   6  * ``Adapted Wavelet Analysis from Theory to Software'' ISBN 1-56881-041-5
   7  * by Mladen Victor Wickerhauser [AK Peters, Ltd., Wellesley, Mass., 1994]
   8  *
   9  */
  10
  11 #include <assert.h>
  12 #include "real.h"
  13 #include "cd.h"
  14 #include "common.h"             /* for `min()' and `max()'. */
  15
  16 /***************************************************************
  17  * cdai():
  18  *
  19  * [C]onvolution and [D]ecimation by 2;
  20  *   [A]periodic, sequential [I]nput, superposition.
  21  *
  22  *  Basic algorithm:
  23  *
  24  *   For J = A to B
  25  *    Let BEGIN = ICH(J+F.ALPHA)
  26  *    Let END   = IFH(J+F.OMEGA)
  27  *    For I = BEGIN to END
  28  *       OUT[I*STEP] += F.F[2*I-J]*IN[J]
  29  *
  30  *  Calling sequence:
  31  *     cdai(out, step, in, a, b, F)
  32  *
  33  *  Inputs:
  34  *    (real *)out   This output array must be preallocated so
  35  *                    that `out[step*ICH(a+F->alpha)]' through
  36  *                    `out[step*IFH(b+F->omega)]' are defined.
  37  *    (int)step     This integer is the increment between
  38  *                    `out[]' elements.
  39  *    (const real *)in    This input array must be preallocated
  40  *                    so that `in[a],...,in[b]' are all defined.
  41  *    (int)a,b      These are the first and last valid indices
  42  *                    for array `in[]'.
  43  *    (const pqf *)F     This specifies the QF to use.
  44  *
  45  *  Outputs:
  46  *    (real *)out   Computed values are added into this array
  47  *                    by side effect, at elements separated
  48  *                    by offsets of `step' starting at
  49  *                    `out + step*ICH(a + F->alpha)'.
  50  */
  51 extern void
  52   cdai(
  53        real *out,      /* Preallocated output array     */
  54        int  step,      /* Increment for `out[]'         */
  55        const real *in, /* Preallocated input array      */
  56        int a,          /* First valid `in[]' offset     */
  57        int b,          /* Last valid `in[]' offset      */
  58        const pqf *F )  /* QF data structure             */
  59 {
  60   int i, j;
  61   real *outp;
  62   const real *filt;
  63
  64   for( j=a; j<=b; j++ )
  65     {
  66       i=ICH(j+F->alpha);
  67       filt = F->f +(2*i-j);
  68       outp = out + i*step;
  69       while( i<= IFH(j+F->omega) )
  70         {
  71           *outp += (*filt) * in[j];
  72           outp += step;   filt += 2;   ++i;
  73         }
  74     }
  75   return;
  76 }
  77 \f
  78 /**************************************************************
  79  * cdao():
  80  *
  81  * [C]onvolution and [D]ecimation by 2;
  82  *   [A]periodic, sequential [O]utput, superposition.
  83  *
  84  *  Basic algorithm:
  85  *
  86  *   Let APRIME = ICH(A+F.ALPHA)
  87  *   Let BPRIME = IFH(B+F.OMEGA)
  88  *   For I = APRIME to BPRIME
  89  *      Let BEGIN = max( F.ALPHA, 2*I-B )
  90  *      Let END   = min( F.OMEGA, 2*I-A )
  91  *      For J = BEGIN to END
  92  *         OUT[I*STEP] += F.F[J]*IN[2*I-J]
  93  *
  94  *  Calling sequence:
  95  *     cdao(out, step, in, a, b, F)
  96  *
  97  *  Inputs:
  98  *    (real *)out   This output array must be preallocated so
  99  *                    that `out[step*ICH(a+F->alpha)]' through
 100  *                    `out[step*IFH(b+F->omega)]' are defined.
 101  *    (int)step     This integer is the increment between
 102  *                    `out[]' elements.
 103  *    (const real *)in    This input array must be preallocated
 104  *                    so that `in[a],...,in[b]' are all defined.
 105  *    (int)a,b      These are the first and last valid indices
 106  *                    for array `in[]'.
 107  *    (const pqf *)F     This specifies the QF to use.
 108  *
 109  *  Outputs:
 110  *    (real *)out   Computed values are added into this array
 111  *                    by side effect, at elements separated
 112  *                    by offsets of `step' starting at
 113  *                    `out + step*ICH(a+F->alpha)'.
 114  */
 115 extern void
 116   cdao(
 117        real *out,      /* Preallocated output array     */
 118        int  step,      /* Increment for `out[]'         */
 119        const real *in, /* Preallocated input array      */
 120        int a,          /* First valid `in[]' offset     */
 121        int b,          /* Last valid `in[]' offset      */
 122        const pqf *F )  /* QF data structure             */
 123 {
 124   int  i, j, end;
 125   const real *inp;
 126
 127   i = ICH(a+F->alpha);
 128   out += i*step;
 129   while( i <= IFH(b+F->omega) )
 130     {
 131       j = 2*i-b;
 132       if(j<F->alpha)
 133         {
 134           j=F->alpha;           /* j = max(F->alpha,2*i-b) */
 135           inp = in + 2*i-j;     /* *inp = in[2*i-j] */
 136         }
 137       else
 138         {
 139           inp = in + b;         /* *inp = in[2*i-j] */
 140         }
 141       end = 2*i-a;
 142       end =  min( F->omega, end );
 143       while( j <= end  )
 144         *out += F->f[j++]*(*inp--);
 145       ++i;
 146       out += step;
 147     }
 148   return;
 149 }
 150 \f
 151 /**************************************************************
 152  * cdae():
 153  *
 154  * [C]onvolution and [D]ecimation by 2;
 155  *   [A]periodic, set output using [E]quals.
 156  *
 157  *  Basic algorithm:
 158  *
 159  *   Let APRIME = ICH(A+F.ALPHA)
 160  *   Let BPRIME = IFH(B+F.OMEGA)
 161  *   For I = APRIME to BPRIME
 162  *      Let OUT[I*STEP] = 0
 163  *      Let BEGIN = max( F.ALPHA, 2*I-B )
 164  *      Let END   = min( F.OMEGA, 2*I-A )
 165  *      For J = BEGIN to END
 166  *         OUT[I*STEP] += F.F[J]*IN[2*I-J]
 167  *
 168  *  Calling sequence:
 169  *     cdae(out, step, in, a, b, F)
 170  *
 171  *  Inputs:
 172  *    (real *)out   This output array must be preallocated so
 173  *                    that `out[step*ICH(a+F->alpha)]' through
 174  *                    `out[step*IFH(b+F->omega)]' are defined.
 175  *    (int)step     This integer is the increment between
 176  *                    `out[]' elements.
 177  *    (const real *)in    This input array must be preallocated
 178  *                    so that `in[a],...,in[b]' are all defined.
 179  *    (int)a,b      These are the first and last valid indices
 180  *                    for array `in[]'.
 181  *    (const pqf *)F     This specifies the QF to use.
 182  *
 183  *  Outputs:
 184  *    (real *)out   Computed values are assigned into this array
 185  *                    by side effect, at elements separated
 186  *                    by offsets of `step' starting at
 187  *                    `out + step*ICH(a+F->alpha)'.
 188  */
 189 extern void
 190   cdae(
 191        real *out,      /* Preallocated output array     */
 192        int  step,      /* Increment for `out[]'         */
 193        const real *in, /* Preallocated input array      */
 194        int a,          /* First valid `in[]' offset     */
 195        int b,          /* Last valid `in[]' offset      */
 196        const pqf *F )  /* QF data structure             */
 197 {
 198   int  i, j, end;
 199   const real *inp;
 200
 201   i = ICH(a+F->alpha);
 202   out += i*step;
 203   while( i <= IFH(b+F->omega) )
 204     {
 205       *out = 0;
 206       j = 2*i-b;
 207       if(j<F->alpha)
 208         {
 209           j=F->alpha;           /* j = max(F->alpha,2*i-b) */
 210           inp = in + 2*i-j;     /* *inp = in[2*i-j] */
 211         }
 212       else
 213         {
 214           inp = in + b;         /* *inp = in[2*i-j] */
 215         }
 216       end = 2*i-a;
 217       end =  min( F->omega, end );
 218       while( j <= end  )
 219         *out += F->f[j++]*(*inp--);
 220       ++i;
 221       out += step;
 222     }
 223   return;
 224 }
 225 \f
 226 /*************************************************************
 227  * acdai():
 228  *
 229  *     [A]djoint [C]onvolution and [D]ecimation:
 230  *      [A]periodic, sequential [I]nput, superposition.
 231  *
 232  *  Basic algorithm:
 233  *
 234  *   For I = C to D
 235  *      Let BEGIN = 2*I-F.OMEGA
 236  *      Let END   = 2*I-F.ALPHA
 237  *      For J = BEGIN to END
 238  *         OUT[J*STEP] += F.F[2*I-J]*IN[I]
 239  *
 240  * We change variables to save operations:
 241  *
 242  *   For I = C to D
 243  *      For J = F.ALPHA to F.OMEGA
 244  *         OUT[(2*I-J)*STEP] += F.F[J]*IN[I]
 245  *
 246  *  Calling sequence:
 247  *     acdai(out, step, in, c, d, F)
 248  *
 249  *  Inputs:
 250  *    (real *)out   This output array must be preallocated so
 251  *                    that `out[step*(2*c-F->omega)]' through
 252  *                    `out[step*(2*d-F->alpha)]' are defined.
 253  *    (int)step     This integer is the increment between
 254  *                    `out[]' elements.
 255  *    (const real *)in    This input array must be preallocated
 256  *                    so that `in[c],...,in[d]' are all defined.
 257  *    (int)a,b      These are the first and last valid indices
 258  *                    for array `in[]'.
 259  *    (const pqf *)F     This specifies the QF to use.
 260  *
 261  *  Outputs:
 262  *    (real *)out   Computed values are added into this array
 263  *                    by side effect, at elements separated
 264  *                    by offsets of `step' starting at
 265  *                    `out + (2*c-F->omega)*step'.
 266  */
 267 extern void
 268   acdai(
 269         real *out,      /* Preallocated output array      */
 270         int  step,      /* Increment for `out[]' values   */
 271         const real *in, /* Preallocated input array       */
 272         int c,          /* First valid `in[]' offset      */
 273         int d,          /* Last valid `in[]' offset       */
 274         const pqf *F )  /* QF data structure              */
 275 {
 276   int i, j;
 277   real *outp;
 278
 279   for( i=c; i <= d ; i++ )
 280     {
 281       j = F->alpha;
 282       outp = out + (2*i-j)*step;
 283       while( j <= F->omega )
 284         {
 285           *outp +=  F->f[j++] * in[i];
 286           outp -= step;
 287         }
 288     }
 289   return;
 290 }
 291 \f
 292 /***************************************************************
 293  * acdao():
 294  *
 295  *     [A]djoint [C]onvolution and [D]ecimation:
 296  *      [A]periodic, sequential [O]utput, superposition.
 297  *
 298  *  Basic algorithm:
 299  *
 300  *   Let CPRIME = 2*C-F.OMEGA
 301  *   Let DPRIME = 2*D-F.ALPHA
 302  *   For J = CPRIME to DPRIME
 303  *      Let BEGIN = max( ICH(J+F.ALPHA ), C )
 304  *      Let END   = min( IFH(J+F.OMEGA ),  D )
 305  *      For I = BEGIN to END
 306  *         OUT[J*STEP] += F.F[2*I-J]*IN[I]
 307  *
 308  *  Calling sequence:
 309  *     acdao(out, step, in, c, d, F)
 310  *
 311  *  Inputs:
 312  *    (real *)out   This output array must be preallocated so
 313  *                    that `out[step*(2*c-F->omega)]' through
 314  *                    `out[step*(2*d-F->alpha)]' are defined.
 315  *    (int)step     This integer is the increment between
 316  *                    `out[]' elements.
 317  *    (const real *)in    This input array must be preallocated
 318  *                    so that `in[c],...,in[d]' are all defined.
 319  *    (int)a,b      These are the first and last valid indices
 320  *                    for array `in[]'.
 321  *    (const pqf *)F     This specifies the QF to use.
 322  *
 323  *  Outputs:
 324  *    (real *)out   Computed values are added into this array
 325  *                    by side effect, at elements separated
 326  *                    by offsets of `step' starting at
 327  *                    `out + (2*c - F->omega)*step'.
 328  */
 329 extern void
 330   acdao(
 331         real *out,      /* Preallocated output array       */
 332         int  step,      /* Increment for `out[]' values    */
 333         const real *in, /* Preallocated input array        */
 334         int c,          /* First valid `in[]' offset       */
 335         int d,          /* Last valid `in[]' offset        */
 336         const pqf *F )  /* QF data structure               */
 337 {
 338   int  end, i, j;
 339   const real *filt;
 340
 341   j = 2*c - F->omega;
 342   out += j*step;
 343
 344   while( j <= 2*d-F->alpha )
 345     {
 346       i = ICH( j + F->alpha );
 347       i = max( i, c );
 348       end   = IFH( j + F->omega );
 349       end   = min( end, d );
 350       filt = F->f+2*i-j;
 351       while( i<=end )
 352         {
 353           *out += (*filt)*in[i++];
 354           filt += 2;
 355         }
 356       out += step;
 357       ++j;
 358     }
 359   return;
 360 }
 361 \f
 362 /***************************************************************
 363  * acdae():
 364  *
 365  *     [A]djoint [C]onvolution and [D]ecimation:
 366  *      [A]periodic, set output using [E]quals.
 367  *
 368  *  Basic algorithm:
 369  *
 370  *   Let CPRIME = 2*C-F.OMEGA
 371  *   Let DPRIME = 2*D-F.ALPHA
 372  *   For J = CPRIME to DPRIME
 373  *      Let OUT[J*STEP] = 0
 374  *      Let BEGIN = max( ICH(J+F.ALPHA ), C )
 375  *      Let END   = min( IFH(J+F.OMEGA ),  D )
 376  *      For I = BEGIN to END
 377  *         OUT[J*STEP] += F.F[2*I-J]*IN[I]
 378  *
 379  *  Calling sequence:
 380  *     acdae(out, step, in, c, d, F)
 381  *
 382  *  Inputs:
 383  *    (real *)out   This output array must be preallocated so
 384  *                    that `out[step*(2*c-F->omega)]' through
 385  *                    `out[step*(2*d-F->alpha)]' are defined.
 386  *    (int)step     This integer is the increment between
 387  *                    `out[]' elements.
 388  *    (const real *)in    This input array must be preallocated
 389  *                    so that `in[c],...,in[d]' are all defined.
 390  *    (int)a,b      These are the first and last valid indices
 391  *                    for array `in[]'.
 392  *    (const pqf *)F     This specifies the QF to use.
 393  *
 394  *  Outputs:
 395  *    (real *)out   Computed values are added into this array
 396  *                    by side effect, at elements separated
 397  *                    by offsets of `step' starting at
 398  *                    `out + (2*c - F->omega)*step'.
 399  */
 400 extern void
 401   acdae(
 402         real *out,      /* Preallocated output array       */
 403         int  step,      /* Increment for `out[]' values    */
 404         const real *in, /* Preallocated input array        */
 405         int c,          /* First valid `in[]' offset       */
 406         int d,          /* Last valid `in[]' offset        */
 407         const pqf *F )  /* QF data structure               */
 408 {
 409   int  end, i, j;
 410   const real *filt;
 411
 412   j = 2*c - F->omega;
 413   out += j*step;
 414
 415   while( j <= 2*d-F->alpha )
 416     {
 417       *out = 0;
 418       i = ICH( j + F->alpha );
 419       i = max( i, c );
 420       end   = IFH( j + F->omega );
 421       end   = min( end, d );
 422       filt = F->f+2*i-j;
 423       while( i<=end )
 424         {
 425           *out += (*filt)*in[i++];
 426           filt += 2;
 427         }
 428       out += step;
 429       ++j;
 430     }
 431   return;
 432 }
 433 \f
 434 /*******************************************************
 435  * cdmi():
 436  *
 437  *   [C]onvolution and [D]ecimation by 2, using the
 438  *      [M]od operator, sequential [I]nput, superposition.
 439  *
 440  *  Basic algorithm:
 441  *
 442  *   Let Q2  = Q/2
 443  *   If Q > F.OMEGA-F.ALPHA then
 444  *      Let FILTER = F.F
 445  *      For J = 0 to Q-1
 446  *         Let JA2 = ICH(J+F.ALPHA)
 447  *         Let JO2 = IFH(J+F.OMEGA)
 448  *         For I = 0 to JO2-Q2
 449  *            OUT[I*STEP] += FILTER[Q+2*I-J]*IN[J]
 450  *         For I = max(0,JA2) to min(Q2-1,JO2)
 451  *            OUT[I*STEP] += FILTER[2*I-J]*IN[J]
 452  *         For I = JA2+Q2 to Q2-1
 453  *            OUT[I*STEP] += FILTER[2*I-J-Q]*IN[J]
 454  *   Else
 455  *      Let FILTER = F.FP+PQFO(Q2)
 456  *      For J = 0 to Q-1
 457  *         For I = 0 to Q2-1
 458  *            OUT[I*STEP] += FILTER[(Q+2*I-J)%Q]*IN[J]
 459  *
 460  *  Calling sequence:
 461  *     cdmi(out, step, in, q, F)
 462  *
 463  *  Inputs:
 464  *    (real *)out   This output array must be preallocated so
 465  *                    that elements `out[0]' through
 466  *                    `out[step*(q/2-1)]' are defined.
 467  *    (int)step     This integer is the increment between
 468  *                    `out[]' elements.
 469  *    (const real *)in    This input array must be preallocated
 470  *                    so that `in[0],...,in[q-1]' are defined.
 471  *    (int)q        This is the period of the input array.
 472  *    (const pqf *)F    This specifies the periodized QF struct.
 473  *
 474  *  Outputs:
 475  *    (real *)out   Computed values are added into this array
 476  *                    by side effect, at elements separated
 477  *                    by offsets of `step' starting at `out'.
 478  *
 479  *  Assumptions: We assume that `q' is even.
 480  */
 481 extern void
 482   cdmi(
 483        real *out,      /* Preallocated array, length `q/2'  */
 484        int step,       /* Increment between `out[]' values  */
 485        const real *in, /* Preallocated array, length `q'    */
 486        int   q,        /* Number of inputs---must be even   */
 487        const pqf *F )  /* Periodized QF data structure      */
 488 {
 489   int i, j, q2;
 490   const real *filt;
 491   real *outp;
 492
 493   assert( (q&1)==0 );           /* Assume `q' is even.  */
 494
 495   q2 = q/2;                     /* Number of outputs.   */
 496   if( q > F->omega-F->alpha )   /* Long input case.     */
 497     {
 498       int  ja2, jo2;
 499
 500       for( j=0; j<q; j++ )
 501         {
 502           ja2 = ICH(j+F->alpha);
 503           jo2 = IFH(j+F->omega);
 504
 505           i=0;    outp = out;    filt = F->f + q-j;
 506           while( i <= jo2-q2 )
 507             {
 508               *outp += (*filt) * (*in);
 509               ++i;  outp += step; filt += 2;
 510             }
 511
 512           i = max(0, ja2);   outp = out+i*step;  filt = F->f + 2*i-j;
 513           while( i <= min(q2-1,jo2) )
 514             {
 515               *outp += (*filt) * (*in);
 516               ++i; outp += step; filt += 2;
 517             }
 518
 519           i = ja2+q2;   outp = out+i*step;   filt = F->f + 2*i-j-q;
 520           while( i < q2 )
 521             {
 522               *outp +=(*filt) * (*in);
 523               ++i; outp += step; filt += 2;
 524             }
 525
 526           ++in;
 527         }
 528     }
 529   else                         /* Short input case.     */
 530     {
 531       filt = F->fp + PQFO(q2);
 532       for( j=0; j<q; j++ )
 533         {
 534           outp = out;
 535           for( i=0; i<q2; i++ )
 536             {
 537               *outp += filt[(q+2*i-j)%q]*(*in);
 538               outp += step;
 539             }
 540           ++in;
 541         }
 542     }
 543   return;
 544 }
 545 \f
 546 /************************************************************
 547  * cdmo():
 548  *
 549  *    [C]onvolution and [D]ecimation by 2, using the
 550  *       [M]od operator,  sequential [O]utput, superposition.
 551  *
 552  *  Basic algorithm:
 553  *
 554  *   Let Q2 = Q/2
 555  *   If Q > F.OMEGA-F.ALPHA then
 556  *      Let FILTER = F.F
 557  *      For I = 0 to Q2-1
 558  *         For J = F.ALPHA to F.OMEGA
 559  *            OUT[I*STEP] += FILTER[J]*IN[(Q+2*I-J)%Q]
 560  *   Else
 561  *      Let FILTER = F.FP + PQFO[Q2]
 562  *      For I = 0 to Q2-1
 563  *         For J = 0 to Q-1
 564  *            OUT[I*STEP] += FILTER[J]*IN[(Q+2*I-J)%Q]
 565  *
 566  *  Calling sequence:
 567  *     cdmo(out, step, in, q, F)
 568  *
 569  *  Inputs:
 570  *    (real *)out   This output array must be preallocated so
 571  *                    that elements `out[0]' through
 572  *                    `out[step*(q/2-1)]' are defined.
 573  *    (int)step     This integer is the increment between
 574  *                    `out[]' elements.
 575  *    (const real *)in    This input array must be preallocated
 576  *                    so that `in[0],...,in[q-1]' are defined.
 577  *    (int)q        This is the period of the input array.
 578  *    (const pqf *)F    This specifies the periodized QF struct.
 579  *
 580  *  Outputs:
 581  *    (real *)out   Computed values are added into this array
 582  *                    by side effect, at elements separated
 583  *                    by offsets of `step' starting at `out'.
 584  *
 585  *  Assumptions: We assume that `q' is even.
 586  */
 587 extern void
 588   cdmo(
 589        real *out,      /* Preallocated array, length `q/2'  */
 590        int step,       /* Increment between `out[]' values  */
 591        const real *in, /* Preallocated array, length `q'    */
 592        int   q,        /* Number of inputs---must be even   */
 593        const pqf *F )  /* Periodized QF data structure      */
 594 {
 595   int i, j, k, q2;
 596
 597   assert( (q&1)==0 );            /* Test that `q' is even.  */
 598
 599   q2 = q/2;
 600   if( q > F->omega - F->alpha )    /* Long input case.    */
 601     {
 602       for( i=0; i<q2; i++)
 603         {
 604           j = F->alpha;
 605           k = q+2*i-j;
 606           while( j <= F->omega)
 607             *out += F->f[j++]*in[(k--)%q];
 608           out += step;
 609         }
 610       }
 611   else                             /* Short input case.   */
 612     {
 613       for(i=0; i<q2; i++)
 614         {
 615           j = PQFO(q2); /* q-periodized coefs. */
 616           k = q+2*i;
 617           while( k > 2*i)
 618             *out += F->fp[j++]*in[(k--)%q];
 619           out += step;
 620         }
 621     }
 622   return;
 623 }
 624 \f
 625 /************************************************************
 626  * cdpi():
 627  *
 628  *    [C]onvolution and [D]ecimation by 2:
 629  *       [P]eriodic, sequential [I]nput, superposition.
 630  *
 631  *  Basic algorithm:
 632  *
 633  *   Let Q2  = Q/2
 634  *   If Q > F.OMEGA-F.ALPHA then
 635  *      Let FILTER = F.F
 636  *      For J = 0 to Q-1
 637  *         Let JA2 = ICH(J+F.ALPHA)
 638  *         Let JO2 = IFH(J+F.OMEGA)
 639  *         For I = 0 to JO2-Q2
 640  *            OUT[I*STEP] += FILTER[Q+2*I-J]*IN[J]
 641  *         For I = max(0,JA2) to min(Q2-1,JO2)
 642  *            OUT[I*STEP] += FILTER[2*I-J]*IN[J]
 643  *         For I = JA2+Q2 to Q2-1
 644  *            OUT[I*STEP] += FILTER[2*I-J-Q]*IN[J]
 645  *   Else
 646  *      Let FILTER = F.FP+PQFO(Q2)
 647  *      For J = 0 to Q-1
 648  *         For I = 0 to Q2-1
 649  *            OUT[I*STEP] += FILTER[2*I-J]*IN[J]
 650  *
 651  *  Calling sequence:
 652  *     cdpi(out, step, in, q, F)
 653  *
 654  *  Inputs:
 655  *    (real *)out   This output array must be preallocated so
 656  *                    that elements `out[0]' through
 657  *                    `out[step*(q/2-1)]' are defined.
 658  *    (int)step     This integer is the increment between
 659  *                    `out[]' elements.
 660  *    (const real *)in    This input array must be preallocated
 661  *                    so that `in[0],...,in[q-1]' are defined.
 662  *    (int)q        This is the period of the input array.
 663  *    (const pqf *)F    This specifies the periodized QF struct.
 664  *
 665  *  Outputs:
 666  *    (real *)out   Computed values are added into this array
 667  *                    by side effect, at elements separated
 668  *                    by offsets of `step' starting at `out'.
 669  *
 670  *  Assumptions: We assume that `q' is even.
 671  */
 672 extern void
 673   cdpi(
 674        real *out,      /* Preallocated array, length `q/2'  */
 675        int step,       /* Increment between `out[]' values  */
 676        const real *in, /* Preallocated array, length `q'    */
 677        int   q,        /* Number of inputs---must be even   */
 678        const pqf *F )  /* Periodized QF data structure      */
 679 {
 680   int i, j, q2;
 681   const real *filt;
 682   real *outp;
 683
 684   assert( (q&1)==0 );          /* Test that `q' is even.   */
 685
 686   q2 = q/2;                    /* Number of outputs.     */
 687   if( q > F->omega-F->alpha )      /* Long input case. */
 688     {
 689       int  ja2, jo2;
 690
 691       for( j=0; j<q; j++ )
 692         {
 693           ja2 = ICH(j+F->alpha);
 694           jo2 = IFH(j+F->omega);
 695
 696           i=0;    outp = out;    filt = F->f + q-j;
 697           while( i<=jo2-q2 )
 698             {
 699               *outp += (*filt) * (*in);
 700               ++i;  outp += step; filt += 2;
 701             }
 702
 703           i = max(0, ja2); outp = out+i*step; filt = F->f + 2*i-j;
 704           while(  i<=min(q2-1,jo2) )
 705             {
 706               *outp += (*filt) * (*in);
 707               ++i; outp += step; filt += 2;
 708             }
 709
 710           i = ja2+q2; outp = out+i*step; filt = F->f + 2*i-j-q;
 711           while( i<q2 )
 712             {
 713               *outp +=(*filt) * (*in);
 714               ++i; outp += step; filt += 2;
 715             }
 716
 717           ++in;
 718         }
 719     }
 720   else                         /* Short input case.     */
 721     {
 722       for( j=0; j<q; j++ )
 723         {
 724           filt = F->fp + PQFO(q2) - j;
 725           outp = out;
 726           for( i=0; i<q2; i++ )
 727             {
 728               *outp += (*filt) * (*in);
 729               filt += 2;   outp += step;
 730             }
 731           ++in;
 732         }
 733     }
 734   return;
 735 }
 736 \f
 737 /*********************************************************
 738  * cdpo1():
 739  *
 740  *     [C]onvolution and [D]ecimation by 2:
 741  *         [P]eriodic, sequential [O]utput, superposition (1).
 742  *
 743  *   Basic algorithm:
 744  *
 745  *   Let Q2  = Q/2
 746  *   If Q > F.OMEGA-F.ALPHA then
 747  *      Let FILTER = F.F
 748  *      For I = 0 to Q2-1
 749  *         Let A2I = 2*I-F.ALPHA
 750  *         Let O2I = 2*I-F.OMEGA
 751  *         For J = 0 to A2I-Q
 752  *            OUT[I*STEP] += FILTER[2*I-J-Q]*IN[J]
 753  *         For J = max(0,O2I) to min(Q-1,A2I)
 754  *            OUT[I*STEP] += FILTER[2*I-J]*IN[J]
 755  *         For J = O2I+Q to Q-1
 756  *            OUT[I*STEP] += FILTER[Q+2*I-J]*IN[J]
 757  *   Else
 758  *      Let FILTER = F.FP+PQFO(Q2)
 759  *      For I = 0 to Q2-1
 760  *         For J = 0 to Q-1
 761  *            OUT[I*STEP] += FILTER[2*I-J]*IN[J]
 762  *
 763  *
 764  *  Calling sequence:
 765  *     cdpo1(out, step, in, q, F)
 766  *
 767  *  Inputs:
 768  *    (real *)out   This output array must be preallocated so
 769  *                    that elements `out[0]' through
 770  *                    `out[step*(q/2-1)]' are defined.
 771  *    (int)step     This integer is the increment between
 772  *                    `out[]' elements.
 773  *    (const real *)in   This input array must be preallocated
 774  *                    so that `in[0],...,in[q-1]' are defined.
 775  *    (int)q        This is the period of the input array.
 776  *    (const pqf *)F     This specifies the periodized QF struct.
 777  *
 778  *  Outputs:
 779  *    (real *)out   Computed values are added into this array
 780  *                    by side effect, at elements separated
 781  *                    by offsets of `step' starting at `out'.
 782  *
 783  *  Assumptions: We assume that `q' is even.
 784  */
 785 extern void
 786   cdpo1(
 787         real *out,      /* Preallocated array, length `q/2'  */
 788         int step,       /* Increment between `out[]' values  */
 789         const real *in, /* Preallocated array, length `q'    */
 790         int   q,        /* Number of inputs---must be even   */
 791         const pqf *F )  /* Periodized QF data structure      */
 792 {
 793   int ii, j;
 794   const real *filt;
 795   int q2;
 796
 797   assert( (q&1)==0 );          /* Test that `q' is even.   */
 798
 799   q2 = q/2;
 800   if( q > F->omega-F->alpha )  /* Long input case.       */
 801     {
 802       int a2i, o2i;
 803
 804       for( ii=0; ii<q; ii+=2 )
 805         {
 806           a2i = ii - F->alpha;
 807           o2i = ii - F->omega;
 808
 809           filt = F->f + ii-q;
 810           for( j=0; j<=a2i-q; j++ )
 811             *out += (*filt--)* in[j];
 812
 813           j = max(0, o2i);  filt = F->f + ii-j;
 814           for( ; j<=min(q-1,a2i); j++ )
 815             *out += (*filt--)* in[j];
 816
 817           j = o2i + q;      filt = F->f + q+ii-j;
 818           for( ; j<q; j++ )
 819             *out += (*filt--) * in[j];
 820
 821           out += step;
 822         }
 823     }
 824   else                         /* Short input case.     */
 825     {
 826       for( ii=0; ii<q; ii+=2 )
 827         {
 828           filt = F->fp + PQFO(q2) + ii;
 829           for( j=0; j<q; j++ )
 830             *out += (*filt--) * in[j];
 831           out += step;
 832         }
 833     }
 834   return;
 835 }
 836 \f
 837 /*********************************************************
 838  * cdpo2():
 839  *
 840  *     [C]onvolution and [D]ecimation by 2:
 841  *         [P]eriodic, sequential [O]utput, superposition (2).
 842  *
 843  *   Basic algorithm:
 844  *
 845  *   Let Q2 = Q/2
 846  *   If Q > F.OMEGA-F.ALPHA then
 847  *      Let FILTER = F.F
 848  *      For I = 0 to Q2-1
 849  *         Let END   = 2*I-Q
 850  *         For J = F.ALPHA to END
 851  *            OUT[I*STEP] += FILTER[J]*IN[2*I-J-Q]
 852  *         Let BEGIN = max(F.ALPHA, END+1)
 853  *         Let END   = min(F.OMEGA, 2*I)
 854  *         For J = BEGIN to END
 855  *            OUT[I*STEP] += FILTER[J]*IN[2*I-J]
 856  *         Let BEGIN = 2*I+1
 857  *         For J = BEGIN to F.OMEGA
 858  *            OUT[I*STEP] += FILTER[J]*IN[Q+2*I-J]
 859  *   Else
 860  *      Let FILTER = F.FP + PQFO(Q2)
 861  *      For I = 0 to Q2-1
 862  *         For J = 0 to Q-1
 863  *            OUT[I*STEP] += FILTER[2*I-J]*IN[J]
 864  *
 865  *  Calling sequence:
 866  *     cdpo2(out, step, in, q, F)
 867  *
 868  *  Inputs:
 869  *    (real *)out   This output array must be preallocated so
 870  *                    that elements `out[0]' through
 871  *                    `out[step*(q/2-1)]' are defined.
 872  *    (int)step     This integer is the increment between
 873  *                    `out[]' elements.
 874  *    (const real *)in   This input array must be preallocated
 875  *                    so that `in[0],...,in[q-1]' are defined.
 876  *    (int)q        This is the period of the input array.
 877  *    (const pqf *)F     This specifies the periodized QF struct.
 878  *
 879  *  Outputs:
 880  *    (real *)out   Computed values are added into this array
 881  *                    by side effect, at elements separated
 882  *                    by offsets of `step' starting at `out'.
 883  *
 884  *  Assumptions: We assume that `q' is even.
 885  */
 886 extern void
 887   cdpo2(
 888         real *out,      /* Preallocated array, length `q/2'  */
 889         int step,       /* Increment between `out[]' values  */
 890         const real *in, /* Preallocated array, length `q'    */
 891         int   q,        /* Number of inputs---must be even   */
 892         const pqf *F )  /* Periodized QF data structure      */
 893 {
 894   int ii, j;
 895   const real *ptr;
 896   int q2;
 897
 898   assert( (q&1)==0 );          /* Test that `q' is even.   */
 899
 900   q2 = q/2;
 901   if( q > F->omega-F->alpha )  /* Long input case.       */
 902     {
 903       int  end;
 904
 905       for( ii=0; ii<q; ii+=2 )
 906         {
 907           end = ii-q;
 908           ptr = in + end;
 909           for( j=F->alpha; j<=end; j++ )
 910             *out += F->f[j]* (*ptr--);
 911
 912           j   = max(F->alpha, j);
 913           end = min(F->omega, ii);
 914           ptr = in + ii - j;
 915           for( ; j<=end; j++ )
 916             *out += F->f[j]* (*ptr--);
 917
 918           j   = ii + 1;
 919           ptr = in + q - 1;
 920           for( ; j<=F->omega; j++ )
 921             *out += F->f[j] * (*ptr--);
 922
 923           out += step;
 924         }
 925     }
 926   else                         /* Short input case.     */
 927     {
 928       for( ii=0; ii<q; ii+=2 )
 929         {
 930           ptr = F->fp + PQFO(q2) + ii;
 931           for( j=0; j<q; j++ )
 932             *out += (*ptr--) * in[j];
 933           out += step;
 934         }
 935     }
 936   return;
 937 }
 938 \f
 939 /*********************************************************
 940  * cdpe1():
 941  *
 942  *     [C]onvolution and [D]ecimation by 2:
 943  *         [P]eriodic, set output using [E]quals (1).
 944  *
 945  *   Basic algorithm:
 946  *
 947  *   Let Q2  = Q/2
 948  *   If Q > F.OMEGA-F.ALPHA then
 949  *      Let FILTER = F.F
 950  *      For I = 0 to Q2-1
 951  *      Let OUT[I*STEP] = 0
 952  *         Let A2I = 2*I-F.ALPHA
 953  *         Let O2I = 2*I-F.OMEGA
 954  *         For J = 0 to A2I-Q
 955  *            OUT[I*STEP] += FILTER[2*I-J-Q]*IN[J]
 956  *         For J = max(0,O2I) to min(Q-1,A2I)
 957  *            OUT[I*STEP] += FILTER[2*I-J]*IN[J]
 958  *         For J = O2I+Q to Q-1
 959  *            OUT[I*STEP] += FILTER[Q+2*I-J]*IN[J]
 960  *   Else
 961  *      Let FILTER = F.FP+PQFO(Q2)
 962  *      For I = 0 to Q2-1
 963  *         Let OUT[I*STEP] = 0
 964  *         For J = 0 to Q-1
 965  *            OUT[I*STEP] += FILTER[2*I-J]*IN[J]
 966  *
 967  *
 968  *  Calling sequence:
 969  *     cdpe1(out, step, in, q, F)
 970  *
 971  *  Inputs:
 972  *    (real *)out   This output array must be preallocated so
 973  *                    that elements `out[0]' through
 974  *                    `out[step*(q/2-1)]' are defined.
 975  *    (int)step     This integer is the increment between
 976  *                    `out[]' elements.
 977  *    (const real *)in   This input array must be preallocated
 978  *                    so that `in[0],...,in[q-1]' are defined.
 979  *    (int)q        This is the period of the input array.
 980  *    (const pqf *)F     This specifies the periodized QF struct.
 981  *
 982  *  Outputs:
 983  *    (real *)out   Computed values are assigned into this array
 984  *                    by side effect, at elements separated
 985  *                    by offsets of `step' starting at `out'.
 986  *
 987  *  Assumptions: We assume that `q' is even.
 988  */
 989 extern void
 990   cdpe1(
 991         real *out,      /* Preallocated array, length `q/2'  */
 992         int step,       /* Increment between `out[]' values  */
 993         const real *in, /* Preallocated array, length `q'    */
 994         int   q,        /* Number of inputs---must be even   */
 995         const pqf *F )  /* Periodized QF data structure      */
 996 {
 997   int ii, j;
 998   const real *filt;
 999   int q2;
1000
1001   assert( (q&1)==0 );          /* Test that `q' is even.   */
1002
1003   q2 = q/2;
1004   if( q > F->omega-F->alpha )  /* Long input case.       */
1005     {
1006       int a2i, o2i;
1007
1008       for( ii=0; ii<q; ii+=2 )
1009         {
1010           *out = 0;
1011           a2i = ii - F->alpha;
1012           o2i = ii - F->omega;
1013
1014           filt = F->f + ii-q;
1015           for( j=0; j<=a2i-q; j++ )
1016             *out += (*filt--)* in[j];
1017
1018           j = max(0, o2i);  filt = F->f + ii-j;
1019           for( ; j<=min(q-1,a2i); j++ )
1020             *out += (*filt--)* in[j];
1021
1022           j = o2i + q;      filt = F->f + q+ii-j;
1023           for( ; j<q; j++ )
1024             *out += (*filt--) * in[j];
1025
1026           out += step;
1027         }
1028     }
1029   else                         /* Short input case.     */
1030     {
1031       for( ii=0; ii<q; ii+=2 )
1032         {
1033           *out = 0;
1034           filt = F->fp + PQFO(q2) + ii;
1035           for( j=0; j<q; j++ )
1036             *out += (*filt--) * in[j];
1037           out += step;
1038         }
1039     }
1040   return;
1041 }
1042 \f
1043 /*********************************************************
1044  * cdpe2():
1045  *
1046  *     [C]onvolution and [D]ecimation by 2:
1047  *         [P]eriodic, output set by [E]qual sign (2).
1048  *
1049  *   Basic algorithm:
1050  *
1051  *   Let Q2 = Q/2
1052  *   If Q > F.OMEGA-F.ALPHA then
1053  *      Let FILTER = F.F
1054  *      For I = 0 to Q2-1
1055  *         Let OUT[I*STEP] = 0
1056  *         Let END   = 2*I-Q
1057  *         For J = F.ALPHA to END
1058  *            OUT[I*STEP] += FILTER[J]*IN[2*I-J-Q]
1059  *         Let BEGIN = max(F.ALPHA, END+1)
1060  *         Let END   = min(F.OMEGA, 2*I)
1061  *         For J = BEGIN to END
1062  *            OUT[I*STEP] += FILTER[J]*IN[2*I-J]
1063  *         Let BEGIN = 2*I+1
1064  *         For J = BEGIN to F.OMEGA
1065  *            OUT[I*STEP] += FILTER[J]*IN[Q+2*I-J]
1066  *   Else
1067  *      Let FILTER = F.FP + PQFO(Q2)
1068  *      For I = 0 to Q2-1
1069  *         Let OUT[I*STEP] = 0
1070  *         For J = 0 to Q-1
1071  *            OUT[I*STEP] += FILTER[2*I-J]*IN[J]
1072  *
1073  *  Calling sequence:
1074  *     cdpe2(out, step, in, q, F)
1075  *
1076  *  Inputs:
1077  *    (real *)out   This output array must be preallocated so
1078  *                    that elements `out[0]' through
1079  *                    `out[step*(q/2-1)]' are defined.
1080  *    (int)step     This integer is the increment between
1081  *                    `out[]' elements.
1082  *    (const real *)in   This input array must be preallocated
1083  *                    so that `in[0],...,in[q-1]' are defined.
1084  *    (int)q        This is the period of the input array.
1085  *    (const pqf *)F     This specifies the periodized QF struct.
1086  *
1087  *  Outputs:
1088  *    (real *)out   Computed values are assigned into this array
1089  *                    by side effect, at elements separated
1090  *                    by offsets of `step' starting at `out'.
1091  *
1092  *  Assumptions: We assume that `q' is even.
1093  */
1094 extern void
1095   cdpe2(
1096         real *out,      /* Preallocated array, length `q/2'  */
1097         int step,       /* Increment between `out[]' values  */
1098         const real *in, /* Preallocated array, length `q'    */
1099         int   q,        /* Number of inputs---must be even   */
1100         const pqf *F )  /* Periodized QF data structure      */
1101 {
1102   int ii, j;
1103   const real *ptr;
1104   int q2;
1105
1106   assert( (q&1)==0 );          /* Test that `q' is even.   */
1107
1108   q2 = q/2;
1109   if( q > F->omega-F->alpha )  /* Long input case.       */
1110     {
1111       int  end;
1112
1113       for( ii=0; ii<q; ii+=2 )
1114         {
1115           *out = 0;
1116           end = ii-q;
1117           ptr = in + end;
1118           for( j=F->alpha; j<=end; j++ )
1119             *out += F->f[j]* (*ptr--);
1120
1121           j   = max(F->alpha, j);
1122           end = min(F->omega, ii);
1123           ptr = in + ii - j;
1124           for( ; j<=end; j++ )
1125             *out += F->f[j]* (*ptr--);
1126
1127           j   = ii + 1;
1128           ptr = in + q - 1;
1129           for( ; j<=F->omega; j++ )
1130             *out += F->f[j] * (*ptr--);
1131
1132           out += step;
1133         }
1134     }
1135   else                         /* Short input case.     */
1136     {
1137       for( ii=0; ii<q; ii+=2 )
1138         {
1139           *out = 0;
1140           ptr = F->fp + PQFO(q2) + ii;
1141           for( j=0; j<q; j++ )
1142             *out += (*ptr--) * in[j];
1143           out += step;
1144         }
1145     }
1146   return;
1147 }
1148 \f
1149 /************************************************************
1150  * acdpi():
1151  *
1152  *    [A]djoint [C]onvolution-[D]ecimation:
1153  *      [P]eriodic, sequential [I]nput, superposition.
1154  *
1155  *  Basic algorithm:
1156  *
1157  *   Let Q = 2*Q2
1158  *   If Q > F.OMEGA-F.ALPHA then
1159  *      Let FILTER = F.F
1160  *      For I = 0 to Q2-1
1161  *         Let A2I = 2*I-F.ALPHA
1162  *         Let O2I = 2*I-F.OMEGA
1163  *         For J = 0 to A2I-Q
1164  *            OUT[J*STEP] += FILTER[2*I-J-Q]*IN[I]
1165  *         For J = max(0,O2I) to min(Q-1,A2I)
1166  *            OUT[J*STEP] += FILTER[2*I-J]*IN[I]
1167  *         For J = O2I+Q to Q-1
1168  *            OUT[J*STEP] += FILTER[Q+2*I-J]*IN[I]
1169  *   Else
1170  *      Let FILTER = F.FP + PQFO(Q2)
1171  *      For I = 0 to Q2-1
1172  *         For J = 0 to Q-1
1173  *            OUT[J*STEP] += FILTER[2*I-J]*IN[I]
1174  *
1175  *
1176  *  Calling sequence:
1177  *     acdpi(out, step, in, q2, F)
1178  *
1179  *  Inputs:
1180  *    (real *)out   This output array must be preallocated so
1181  *                    that elements `out[0]' through
1182  *                    `out[step*(q-1)]' are defined, q=2*q2.
1183  *    (int)step     This integer is the increment between
1184  *                    `out[]' elements.
1185  *    (const real *)in    This input array must be preallocated
1186  *                    so that `in[0],...,in[q2-1]' are defined.
1187  *    (int)q2       This is the period of the input array.
1188  *    (const pqf *)F     This specifies the periodized QF struct.
1189  *
1190  *  Outputs:
1191  *    (real *)out   Computed values are added into this array
1192  *                    by side effect, at elements separated
1193  *                    by offsets of `step' starting at `out'.
1194  */
1195 extern void
1196   acdpi(
1197         real *out,      /* Preallocated array, length `2*q2' */
1198         int step,       /* Increment between `out[]' values  */
1199         const real *in, /* Preallocated array, length `q2'   */
1200         int q2,         /* Number of elements of `in[]'      */
1201         const pqf *F )  /* Periodized QF data structure      */
1202 {
1203   int ii, j, q;
1204   real *outp;
1205   const real *filt;
1206
1207   q = 2*q2;                      /* Output array length.  */
1208   if( q > F->omega-F->alpha )    /* Long signal case:     */
1209     {
1210       int a2i, o2i;
1211
1212       for( ii=0; ii<q; ii+=2 )
1213         {
1214           a2i = ii - F->alpha;
1215           o2i = ii - F->omega;
1216
1217           j = 0;   outp = out;   filt = F->f + ii-q;
1218           while( j<=a2i-q )
1219             {
1220               *outp += (*filt--) * (*in);
1221               outp += step; ++j;
1222             }
1223
1224           j = max(0, o2i);  outp = out+j*step;  filt = F->f +ii-j;
1225           while( j<=min(a2i, q-1) )
1226             {
1227               *outp += (*filt--) * (*in);
1228               outp += step; ++j;
1229             }
1230
1231           j = o2i+q;  outp = out+j*step;  filt = F->f + q+ii-j;
1232           while( j<q )
1233             {
1234               *outp += (*filt--) * (*in);
1235               outp += step; ++j;
1236             }
1237           ++in;
1238         }
1239     }
1240   else                          /* Short signal case:    */
1241     {
1242       for( ii=0; ii<q; ii+=2 )
1243         {
1244           outp = out;   filt = F->fp + PQFO(q2) + ii;
1245           for( j=0; j<q; j++)
1246             {
1247               *outp += (*filt--) * (*in);
1248               outp += step;
1249             }
1250           ++in;
1251         }
1252     }
1253   return;
1254 }
1255 \f
1256 /*************************************************************
1257  * acdpo():
1258  *
1259  *   [A]djoint [C]onvolution-[D]ecimation:
1260  *       [P]eriodic, sequential [O]utput, superposition.
1261  *
1262  *  Basic algorithm:
1263  *
1264  *   Let Q = 2*Q2
1265  *   If Q > F.OMEGA-F.ALPHA then
1266  *      Let FILTER = F.F
1267  *      For J = 0 to Q-1
1268  *         Let JA2 = ICH(J+F.ALPHA)
1269  *         Let JO2 = IFH(J+F.OMEGA)
1270  *         For I = 0 to JO2-Q2
1271  *            OUT[J*STEP] += FILTER[Q+2*I-J]*IN[I]
1272  *         For I = max(0,JA2) to min(Q2-1,JO2)
1273  *            OUT[J*STEP] += FILTER[2*I-J]*IN[I]
1274  *         For I = JA2+Q2 to Q2-1
1275  *            OUT[J*STEP] += FILTER[2*I-J-Q]*IN[I]
1276  *   Else
1277  *      Let FILTER = F.FP+PQFO(Q2)
1278  *      For J = 0 to Q-1
1279  *         For I = 0 to Q2-1
1280  *            OUT[J*STEP] += FILTER[2*I-J]*IN[I]
1281  *
1282  *  Calling sequence:
1283  *     acdpo(out, step, in, q2, F)
1284  *
1285  *  Inputs:
1286  *    (real *)out   This output array must be preallocated so
1287  *                    that elements `out[0]' through
1288  *                    `out[step*(q-1)]' are defined, q=2*q2.
1289  *    (int)step     This integer is the increment between
1290  *                    `out[]' elements.
1291  *    (const real *)in    This input array must be preallocated
1292  *                    so that `in[0],...,in[q2-1]' are defined.
1293  *    (int)q2       This is the period of the input array.
1294  *    (const pqf *)F     This specifies the periodized QF struct.
1295  *
1296  *  Outputs:
1297  *    (real *)out   Computed values are added into this array
1298  *                    by side effect, at elements separated
1299  *                    by offsets of `step' starting at `out'.
1300  */
1301 extern void
1302   acdpo(
1303         real *out,      /* Preallocated array, length `2*q2' */
1304         int step,       /* Increment between `out[]' values  */
1305         const real *in, /* Preallocated array, length `q2'   */
1306         int   q2,       /* Number of elements of `in[]'      */
1307         const pqf *F )  /* Periodized qf data structure      */
1308 {
1309   int i, j, q;
1310   const real *inp;
1311   const real *filt;
1312
1313   q = 2*q2;                      /* Output array length.   */
1314   if( q > F->omega-F->alpha )    /* Long signal case:      */
1315     {
1316       int ja2, jo2;
1317
1318       for( j=0; j<q; j++ )
1319         {
1320           ja2 = ICH(j+F->alpha);
1321           jo2 = IFH(j+F->omega);
1322
1323           i=0;  filt = F->f +q-j;  inp = in;
1324           while( i<=jo2-q2 )
1325             {
1326               *out += (*filt) * (*inp++);
1327               filt +=2;   ++i;
1328             }
1329
1330           i=max(0, ja2);  filt = F->f + 2*i-j; inp = in+i;
1331           while( i<=min(q2-1,jo2) )
1332             {
1333               *out += (*filt) * (*inp++);
1334               filt +=2;   ++i;
1335             }
1336
1337           i=q2+ja2;  filt = F->f+2*i-j-q;  inp = in+i;
1338           while( i<q2 )
1339             {
1340               *out += (*filt) * (*inp++);
1341               filt +=2;   ++i;
1342             }
1343           out += step;
1344         }
1345     }
1346   else                          /* Short signal case:    */
1347     {
1348       for( j=0; j<q; j++ )
1349         {
1350           inp = in;  filt = F->fp + PQFO(q2) - j;
1351           for( i=0; i<q2; i++)
1352             {
1353               *out += (*filt) * (*inp++);
1354               filt += 2;
1355             }
1356           out += step;
1357         }
1358     }
1359   return;
1360 }
1361 \f
1362 /*************************************************************
1363  * acdpe():
1364  *
1365  *   [A]djoint [C]onvolution-[D]ecimation:
1366  *       [P]eriodic, assign output with [E]quals.
1367  *
1368  *  Basic algorithm:
1369  *
1370  *   Let Q = 2*Q2
1371  *   If Q > F.OMEGA-F.ALPHA then
1372  *      Let FILTER = F.F
1373  *      Let OUT[J*STEP] = 0
1374  *      For J = 0 to Q-1
1375  *         Let JA2 = ICH(J+F.ALPHA)
1376  *         Let JO2 = IFH(J+F.OMEGA)
1377  *         For I = 0 to JO2-Q2
1378  *            OUT[J*STEP] += FILTER[Q+2*I-J]*IN[I]
1379  *         For I = max(0,JA2) to min(Q2-1,JO2)
1380  *            OUT[J*STEP] += FILTER[2*I-J]*IN[I]
1381  *         For I = JA2+Q2 to Q2-1
1382  *            OUT[J*STEP] += FILTER[2*I-J-Q]*IN[I]
1383  *   Else
1384  *      Let FILTER = F.FP+PQFO(Q2)
1385  *      For J = 0 to Q-1
1386  *         Let OUT[J*STEP] = 0
1387  *         For I = 0 to Q2-1
1388  *            OUT[J*STEP] += FILTER[2*I-J]*IN[I]
1389  *
1390  *  Calling sequence:
1391  *     acdpe(out, step, in, q2, F)
1392  *
1393  *  Inputs:
1394  *    (real *)out   This output array must be preallocated so
1395  *                    that elements `out[0]' through
1396  *                    `out[step*(q-1)]' are defined, q=2*q2.
1397  *    (int)step     This integer is the increment between
1398  *                    `out[]' elements.
1399  *    (const real *)in    This input array must be preallocated
1400  *                    so that `in[0],...,in[q2-1]' are defined.
1401  *    (int)q2       This is the period of the input array.
1402  *    (const pqf *)F     This specifies the periodized QF struct.
1403  *
1404  *  Outputs:
1405  *    (real *)out   Computed values are assigned into this array
1406  *                    by side effect, at elements separated
1407  *                    by offsets of `step' starting at `out'.
1408  */
1409 extern void
1410   acdpe(
1411         real *out,      /* Preallocated array, length `2*q2' */
1412         int step,       /* Increment between `out[]' values  */
1413         const real *in, /* Preallocated array, length `q2'   */
1414         int   q2,       /* Number of elements of `in[]'      */
1415         const pqf *F )  /* Periodized QF data structure      */
1416 {
1417   int i, j, q;
1418   const real *inp;
1419   const real *filt;
1420
1421   q = 2*q2;                      /* Output array length.   */
1422   if( q > F->omega-F->alpha )    /* Long signal case:      */
1423     {
1424       int ja2, jo2;
1425
1426       for( j=0; j<q; j++ )
1427         {
1428           *out = 0;
1429           ja2 = ICH(j+F->alpha);
1430           jo2 = IFH(j+F->omega);
1431
1432           i=0;  filt = F->f +q-j;  inp = in;
1433           while( i<=jo2-q2 )
1434             {
1435               *out += (*filt) * (*inp++);
1436               filt +=2;   ++i;
1437             }
1438
1439           i=max(0, ja2);  filt = F->f + 2*i-j; inp = in+i;
1440           while( i<=min(q2-1,jo2) )
1441             {
1442               *out += (*filt) * (*inp++);
1443               filt +=2;   ++i;
1444             }
1445
1446           i=q2+ja2;  filt = F->f+2*i-j-q;  inp = in+i;
1447           while( i<q2 )
1448             {
1449               *out += (*filt) * (*inp++);
1450               filt +=2;   ++i;
1451             }
1452           out += step;
1453         }
1454     }
1455   else                          /* Short signal case:    */
1456     {
1457       for( j=0; j<q; j++ )
1458         {
1459           *out = 0;
1460           inp = in;  filt = F->fp + PQFO(q2) - j;
1461           for( i=0; i<q2; i++)
1462             {
1463               *out += (*filt) * (*inp++);
1464               filt += 2;
1465             }
1466           out += step;
1467         }
1468     }
1469   return;
1470 }