## DF01MD

### Sine transform or cosine transform of a real signal

[Specification] [Arguments] [Method] [References] [Comments] [Example]

Purpose

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

```
Specification
```      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(*)

```
Arguments

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.

```
Input/Output Parameters
```  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.

```
Workspace
```  DWORK   DOUBLE PRECISION array, dimension (N+1)

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

```
Method
```  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

```
References
```   Rabiner, L.R. and Rader, C.M.
Digital Signal Processing.
IEEE Press, 1972.

 Oppenheim, A.V. and Schafer, R.W.
Discrete-Time Signal Processing.
Prentice-Hall Signal Processing Series, 1989.

```
Numerical Aspects
```  The algorithm requires 0( N*log(N) ) operations.

```
```  None
```
Example

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
```
Program Data
``` 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
```
Program Results
``` 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
```