**Purpose**

To compute the sine transform or cosine transform of a real signal.

SUBROUTINE DF01MD( SICO, N, DT, A, DWORK, INFO ) C .. Scalar Arguments .. CHARACTER SICO INTEGER INFO, N DOUBLE PRECISION DT C .. Array Arguments .. DOUBLE PRECISION A(*), DWORK(*)

**Mode Parameters**

SICO CHARACTER*1 Indicates whether the sine transform or cosine transform is to be computed as follows: = 'S': The sine transform is computed; = 'C': The cosine transform is computed.

N (input) INTEGER The number of samples. N must be a power of 2 plus 1. N >= 5. DT (input) DOUBLE PRECISION The sampling time of the signal. A (input/output) DOUBLE PRECISION array, dimension (N) On entry, this array must contain the signal to be processed. On exit, this array contains either the sine transform, if SICO = 'S', or the cosine transform, if SICO = 'C', of the given signal.

DWORK DOUBLE PRECISION array, dimension (N+1)

INFO INTEGER = 0: successful exit; < 0: if INFO = -i, the i-th argument had an illegal value.

Let A(1), A(2),..., A(N) be a real signal of N samples. If SICO = 'S', the routine computes the sine transform of A as follows. First, transform A(i), i = 1,2,...,N, into the complex signal B(i), i = 1,2,...,(N+1)/2, where B(1) = -2*A(2), B(i) = {A(2i-2) - A(2i)} - j*A(2i-1) for i = 2,3,...,(N-1)/2, B((N+1)/2) = 2*A(N-1) and j**2 = -1. Next, perform a discrete inverse Fourier transform on B(i) by calling SLICOT Library Routine DG01ND, to give the complex signal Z(i), i = 1,2,...,(N-1)/2, from which the real signal C(i) may be obtained as follows: C(2i-1) = Re(Z(i)), C(2i) = Im(Z(i)) for i = 1,2,...,(N-1)/2. Finally, compute the sine transform coefficients S ,S ,...,S 1 2 N given by S = 0, 1 { [C(k) + C(N+1-k)] } S = DT*{[C(k) - C(N+1-k)] - -----------------------}, k { [2*sin(pi*(k-1)/(N-1))]} for k = 2,3,...,N-1, and S = 0. N If SICO = 'C', the routine computes the cosine transform of A as follows. First, transform A(i), i = 1,2,...,N, into the complex signal B(i), i = 1,2,...,(N+1)/2, where B(1) = 2*A(1), B(i) = 2*A(2i-1) + 2*j*{[A(2i-2) - A(2i)]} for i = 2,3,...,(N-1)/2 and B((N+1)/2) = 2*A(N). Next, perform a discrete inverse Fourier transform on B(i) by calling SLICOT Library Routine DG01ND, to give the complex signal Z(i), i = 1,2,...,(N-1)/2, from which the real signal D(i) may be obtained as follows: D(2i-1) = Re(Z(i)), D(2i) = Im(Z(i)) for i = 1,2,...,(N-1)/2. Finally, compute the cosine transform coefficients S ,S ,...,S 1 2 N given by S = 2*DT*[D(1) + A0], 1 { [D(k) - D(N+1-k)] } S = DT*{[D(k) + D(N+1-k)] - -----------------------}, k { [2*sin(pi*(k-1)/(N-1))]} for k = 2,3,...,N-1, and S = 2*DT*[D(1) - A0], N (N-1)/2 where A0 = 2*SUM A(2i). i=1

[1] Rabiner, L.R. and Rader, C.M. Digital Signal Processing. IEEE Press, 1972. [2] Oppenheim, A.V. and Schafer, R.W. Discrete-Time Signal Processing. Prentice-Hall Signal Processing Series, 1989.

The algorithm requires 0( N*log(N) ) operations.

None

**Program Text**

* DF01MD EXAMPLE PROGRAM TEXT * Copyright (c) 2002-2017 NICONET e.V. * * .. Parameters .. INTEGER NIN, NOUT PARAMETER ( NIN = 5, NOUT = 6 ) INTEGER NMAX PARAMETER ( NMAX = 129 ) * .. Local Scalars .. DOUBLE PRECISION DT INTEGER I, INFO, N CHARACTER*1 SICO * .. Local Arrays .. DOUBLE PRECISION A(NMAX), DWORK(NMAX+1) * .. External Functions .. LOGICAL LSAME EXTERNAL LSAME * .. External Subroutines .. EXTERNAL DF01MD * .. Executable Statements .. * WRITE ( NOUT, FMT = 99999 ) * Skip the heading in the data file and read the data. READ ( NIN, FMT = '()' ) READ ( NIN, FMT = * ) N, DT, SICO IF ( N.LE.1 .OR. N.GT.NMAX ) THEN WRITE ( NOUT, FMT = 99994 ) N ELSE READ ( NIN, FMT = * ) ( A(I), I = 1,N ) * Compute the sine/cosine transform of the given real signal. CALL DF01MD( SICO, N, DT, A, DWORK, INFO ) * IF ( INFO.NE.0 ) THEN WRITE ( NOUT, FMT = 99998 ) INFO ELSE IF ( LSAME( SICO, 'S' ) ) THEN WRITE ( NOUT, FMT = 99997 ) DO 20 I = 1, N WRITE ( NOUT, FMT = 99995 ) I, A(I) 20 CONTINUE ELSE WRITE ( NOUT, FMT = 99996 ) DO 40 I = 1, N WRITE ( NOUT, FMT = 99995 ) I, A(I) 40 CONTINUE END IF END IF END IF * STOP * 99999 FORMAT (' DF01MD EXAMPLE PROGRAM RESULTS',/1X) 99998 FORMAT (' INFO on exit from DF01MD = ',I2) 99997 FORMAT (' Components of sine transform are',//' i',6X,'A(i)',/) 99996 FORMAT (' Components of cosine transform are',//' i',6X,'A(i)', $ /) 99995 FORMAT (I4,3X,F8.4) 99994 FORMAT (/' N is out of range.',/' N = ',I5) END

DF01MD EXAMPLE PROGRAM DATA 17 1.0 C -0.1862 0.1288 0.3948 0.0671 0.6788 -0.2417 0.1861 0.8875 0.7254 0.9380 0.5815 -0.2682 0.4904 0.9312 -0.9599 -0.3116 0.8743

DF01MD EXAMPLE PROGRAM RESULTS Components of cosine transform are i A(i) 1 28.0536 2 3.3726 3 -20.8158 4 6.0566 5 5.7317 6 -3.9347 7 -12.8074 8 -6.8780 9 16.2892 10 -17.0788 11 21.7836 12 -20.8203 13 -7.3277 14 -2.5325 15 -0.3636 16 7.8792 17 11.0048