dffti.src

来自「没有说明」· SRC 代码 · 共 95 行

/*
** dffti.src
** (C) Copyright 1988-1998 by Aptech Systems, Inc.
** All Rights Reserved.
**
** This Software Product is PROPRIETARY SOURCE CODE OF APTECH
** SYSTEMS, INC.    This File Header must accompany all files using
** any portion, in whole or in part, of this Source Code.   In
** addition, the right to create such files is strictly limited by
** Section 2.A. of the GAUSS Applications License Agreement
** accompanying this Software Product.
**
** If you wish to distribute any portion of the proprietary Source
** Code, in whole or in part, you must first obtain written
** permission from Aptech Systems.
**
**> dffti
**
**  Purpose:    Computes an inverse complex discrete Fourier transform.
**
**      ( GAUSS-386i complex version )
**
**  Format:     y = dffti(x);
**
**  Input:      x    Nx1 complex vector.
**
**  Output:     y    Nx1 complex vector.
**
**      ( GAUSS-386 real version )
**
**  Format:     { yr,yi } = dffti(xr,xi);
**
**  Input:      xr    Nx1 vector, real part.
**
**              xi    Nx1 vector, imaginary part.
**
**  Output:     yr    Nx1 vector, real part.
**
**              yi    Nx1 vector, imaginary part.
**
**  Remarks:    The transform is divided by N.
**
**               This uses a second-order Goertzel algorithm. It is
**               considerably slower than either ffti or rffti, but
**               it may have some advantages in some
**               circumstances. For one thing, N does not have to
**               be an even power of 2.
**
**  Globals:    None
**
**  See Also:   fft, rfft, dfft, ffti, rffti
*/

#ifcplx

proc (1) = dffti(x);
    local N,xy,q,ab,k,s,c,cc1,z2,qk,e2,z;
    z = real(x) ~ imag(x);

#else

proc (2) = dffti(xr,xi);
    local N,xy,q,ab,k,s,c,cc1,z2,qk,e2,z;
    z = xr~xi;

#endif

    N = rows(z);
    q = 2*pi/N;
    ab = zeros(N,2);
    z2 = zeros(2,2);
    e2 = -eye(2);
    xy = z2|z;
    k = 1;
    do until k > N;
        qk = q*(1-k);
        s = sin(qk);
        c = cos(qk);
        cc1 = reshape( ( (2*c)|-1 ), 2, 2)';
        ab[k,.] = reshape(trimr(recserar(xy,z2,cc1),n,0),4,1)'(e2|c~s|-s~c);
        k = k + 1;
    endo;

#ifcplx

    retp( complex(ab[.,1],ab[.,2]) );

#else

    retp(ab[.,1],ab[.,2]);

#endif

endp;