CFIRM1D(3S)CFIRM1D(3S)NAME
CFIRM1D, ZFIRM1D, SFIRM1D, DFIRM1D - Compute multiple 1D convolutions
SYNOPSIS
Single precision complex
Fortran:
CALL CFIRM1D (x, incx, ldx, ix0, nx, nseq, h, inch, ih0, nh, y,
incy, ldy, iy0, ny, alpha, beta)
C/C++:
#include <scsl_fft.h>
void cfirm1d( scsl_complex *x, int incx, int ix0, int nx,
scsl_complex *h, int inch, int ih0, int nh, scsl_complex *y,
int incy, int iy0, int ny, scsl_complex *alpha, scsl_complex
*beta);
C++ STL:
#include <complex.h>
#include <scsl_fft.h>
void cfirm1d( complex<float> *x, int incx, int ix0, int nx,
complex<float> *h, int inch, int ih0, int nh, complex<float>
*y, int incy, int iy0, int ny, complex<float> *alpha,
complex<float> *beta)
Double precision complex
Fortran:
CALL ZFIRM1D (x, incx, ldx, ix0, nx, nseq, h, inch, ih0, nh, y,
incy, ldy, iy0, ny, alpha, beta)
C/C++:
#include <scsl_fft.h>
void zfirm1d( scsl_zomplex *x, int incx, int ix0, int nx,
scsl_zomplex *h, int inch, int ih0, int nh, scsl_zomplex *y,
int incy, int iy0, int ny, scsl_zomplex *alpha, scsl_zomplex
*beta);
C++ STL:
#include <complex.h>
#include <scsl_fft.h>
void zfirm1d( complex<double> *x, int incx, int ix0, int nx,
complex<double> *h, int inch, int ih0, int nh, complex<double>
*y, int incy, int iy0, int ny, complex<double> *alpha,
complex<double> *beta);
Single precision
Fortran:
CALL SFIRM1D (x, incx, ldx, ix0, nx, nseq, h, inch, ih0, nh, y,
incy, ldy, iy0, ny, alpha, beta)
Page 1
CFIRM1D(3S)CFIRM1D(3S)
C/C++:
#include <scsl_fft.h>
void sfirm1d( float *x, int incx, int ix0, int nx, float *h,
int inch, int ih0, int nh, float *y, int incy, int iy0, int ny,
float alpha, float beta);
Double precision
Fortran:
CALL DFIRM1D (x, incx, ldx, ix0, nx, nseq, h, inch, ih0, nh, y,
incy, ldy, iy0, ny, alpha, beta)
C/C++:
#include <scsl_fft.h>
void dfirm1d( double *x, int incx, int ix0, int nx, double *h,
int inch, int ih0, int nh, double *y, int incy, int iy0, int
ny, double alpha, double beta);
IMPLEMENTATION
These routines are part of the SCSL Scientific Library and can be loaded
using either the -lscs or the -lscs_mp option. The -lscs_mp option
directs the linker to use the multi-processor version of the library.
When linking to SCSL with -lscs or -lscs_mp, the default integer size is
4 bytes (32 bits). Another version of SCSL is available in which integers
are 8 bytes (64 bits). This version allows the user access to larger
memory sizes and helps when porting legacy Cray codes. It can be loaded
by using the -lscs_i8 option or the -lscs_i8_mp option. A program may use
only one of the two versions; 4-byte integer and 8-byte integer library
calls cannot be mixed.
The C and C++ prototypes shown above are appropriate for the 4-byte
integer version of SCSL. When using the 8-byte integer version, the
variables of type int become long long and the <scsl_fft_i8.h> header
file should be included.
DESCRIPTION
These routines compute the convolutions of the filter vector h with each
column of the two-dimenional array x, producing the output two-
dimensional array y:
y = beta * y + alpha * h * x
Suppose h is a sequence of nh elements and X is a 2D matrix with nseq
columns, and nx elements in each column, as follows:
h = [ h(0), h(1), , h(nh - 1) ] ,
and
x(0, 0) x(0, 1) x(0, 2) x(0, nseq-1)
x(1, 0) x(1, 1) x(1, 2) x(1, nseq-1)
Page 2
CFIRM1D(3S)CFIRM1D(3S)
X = x(2, 0) x(2, 1) x(2, 2) x(2, nseq-1)
... ... ... ... ...
x(nx-1, 0) x(nx-1, 1) x(nx-1, 2) x(nx-1, nseq-1)
Then each column of the output matrix:
y(0, 0) y(0, 1) y(0, 2) y(0, nseq-1)
y(1, 0) y(1, 1) y(1, 2) y(1, nseq-1)
Y = y(2, 0) y(2, 1) y(2, 2) y(2, nseq-1)
... ... ... ... ...
y(ny-1, 0) y(ny-1, 1) y(ny-1, 2) y(ny-1, nseq-1)
is obtained by convolving h with the corresponding column of so that:
MIN(i, nh-1)
y(i,j)= Sum {h(k)*x(i-k,j)}
k = MAX(0, i-nx+1)
nh-1
y(i,j)= Sum {h(k)*x(i-k,j)}
k=0
0 <=i<nx, 0<=j<nseq
That is:
y(0,j)=h(0)*x(0,j)
y(1,j)=h(0)*x(1,j) + h(1)*x(0,j)
y(2,j)=h(0)*x(2,j) + h(1)*x(1,j)+h(2)*x(0,j)
y(nh-1,j)=h(0)*x(nh-1,j) + h(1)*x(nh-2,j) + ... + h(nh-1)*x(0,j)
...
y(k,j)=h(0)*x(k,j)+h(1)*x(k,j-1) + ... + h(nh-1)*x(k,j-nh+1)
...
y(nx-1,j)=h(0)*x(nx-1,j) + h(1)*x(nx-2,j) + ... + h(nh-1)*x(nx-nh,j)
...
y(nx+nh-3,j)=h(nh-2)*x(nx-1,j) + h(nh-1)*x(nx-2,j)
y(nx+nh-2,j)=h(nh-1)*x(nx-1,j)
In the *FIRM1D routines, the number of terms in the each output column is
specified by an argument, ny. If ny < nh + nx - 1 the columns of y are
truncated. If ny > nh + nx - 1 the terms beyond y(nh + nx - 2) are set
to 0.
Generally, the sequences x(:,j), h and y(:,j) represent signals sampled
at equal time intervals, and the indexes of the vectors denote the sample
times. If the signals begin at the same time, we may, without loss of
generality, set the initial time to 0, as in the formulas above.
Page 3
CFIRM1D(3S)CFIRM1D(3S)
The *FIRM1D routines, however, permit more generality than this. The
signals may be time shifted from each other using input parameters
specifiying the initial time sample for each signal. This can be useful
in several situations. For example, if the input array has several
leading zero values that one does not wish to store, ix0 may be set to
the time sample corresponding to the first non-zero element in the input
array, and earlier time samples are treated as 0. Another use is to limit
the output to just the "fully engaged" terms of the convolution.
As can be seen above, when nx>= nh, the convolution has ramp-up and
ramp-down regions in which fewer than all nh filter values contribute to
the output value. Setting iy0 to nh-1 causes the first value output to
correspond to time sample nh-1, thus skipping the ramp-up region.
Setting ny to nx-nh+1 then drops the ramp-down terms, limiting the output
to just the fully engaged part.
Note that, instead of 0, the initial time could just as easily have been
labeled 1 or 10 or -78; the relevant point is that the first elements of
each of the x, h and y arrays are defined to be the same time sample as
long as ix0 = ih0 = iy0.
See the NOTES section of this man page for information about the
interpretation of the data types described in the following arguments.
These routines have the following arguments:
x Array of dimension (ldx, nseq). (input).
CFIRM1D: Single precision complex array.
ZFIRM1D: Double precision complex array.
SFIRM1D: Single precision array.
DFIRM1D: Double precision array.
Input sequences to be correlated with h.
incx Integer. (input)
Increment between two successive values of a sequence in x.
incx must not be 0.
ldx Integer. (input)
The number of rows in x as it was declared in the calling
program (the leading dimension of x). ldx >= MAX(nx *
incx,1).
ix0 Integer. (input)
Time sample corresponding to the first element of each 1D
sequence of x.
nx Integer. (input)
The number of elements in each sequence of x. nx >= 0.
nseq Integer. (input)
The number of sequences to which the convolution will be
applied. nseq >= 0. If nseq = 0, the routine returns.
Page 4
CFIRM1D(3S)CFIRM1D(3S)
h Array of dimension nh. (input).
CFIRM1D: Single precision complex array.
ZFIRM1D: Double precision complex array.
SFIRM1D: Single precision array.
DFIRM1D: Double precision array.
Input sequence to be convoluted with x.
inch Integer. (input)
Increment between two successive values of h. inch must not be
0.
ih0 Integer. (input)
Time sample corresponding to the first element of h.
nh Integer. (input)
The number of elements in the sequence h. nh >= 0. If nh = 0,
the routine returns.
y Array dimensioned (ldy, nseq). (output)
CFIRM1D: Single precision complex array.
ZFIRM1D: Double precision complex array.
SFIRM1D: Single precision array.
DFIRM1D: Double precision array.
Output of the FIR filter. On entry, the array y must have been
initialized, except when beta is zero. In that case, y need
not be initialized. On exit, the result overwrites y.
incy Integer. (input)
Increment between two successive values of a sequence in y.
incy must not be 0.
ldy Integer. (input)
The number of rows in y as it was declared in the calling
program (the leading dimension of y). ldy >= MAX(ny * incy, 1).
iy0 Integer. (input)
Time sample corresponding to the first element of each 1D
sequence of y.
ny Integer. (input)
Number of elements in each sequence of y. ny >= 0. If ny = 0,
the routine returns.
alpha Scale factor for the convolution. (input).
CFIRM1D: Single precision complex.
ZFIRM1D: Double precision complex.
SFIRM1D: Single precision.
DFIRM1D: Double precision.
For C/C++, a pointer to this value is passed.
Page 5
CFIRM1D(3S)CFIRM1D(3S)
beta Scale factor for the output y. (input)
CFIRM1D: Single precision complex.
ZFIRM1D: Double precision complex.
SFIRM1D: Single precision.
DFIRM1D: Double precision.
When beta is supplied as 0, y need not be set on input. For
C/C++, a pointer to this value is passed.
NOTES
The following data types are described in this documentation:
Term Used Data type
Fortran:
Array dimensioned 0..n-1 x(0:n-1)
Array of dimensions (m,n) x(m,n)
Array of dimensions (m,n,p) x(m,n,p)
Integer INTEGER (INTEGER*8 for -lscs_i8[_mp])
Single precision REAL
Double precision DOUBLE PRECISION
Single precision complex COMPLEX
Double precision complex DOUBLE COMPLEX
C/C++:
Array dimensioned 0..n-1 x[n]
Array of dimensions (m,n) x[m*n] or x[n][m]
Array of dimensions (m,n,p) x[m*n*p] or x[p][n][m]
Integer int (long long for -lscs_i8[_mp])
Single precision float
Double precision double
Single precision complex scsl_complex
Double precision complex scsl_zomplex
C++ STL:
Page 6
CFIRM1D(3S)CFIRM1D(3S)
Array dimensioned 0..n-1 x[n]
Array of dimensions (m,n) x[m*n] or x[n][m]
Array of dimensions (m,n,p) x[m*n*p] or x[p][n][m]
Integer int (long long for -lscs_i8[_mp])
Single precision float
Double precision double
Single precision complex complex<float>
Double precision complex complex<double>
CAUTIONS
The arrays x, h, and y must be non-overlapping.
EXAMPLES
The following example computes the convolution of 5 4-sample sequences x
with a filter h containing 3 samples:
Fortran:
REAL X(0:3,0:4), H(0:2), Y(0:5,0:4)
REAL ALPHA, BETA
ALPHA = 1.0
BETA = 0.0
DO J = 0, 4
X(0,J) = J + 1.0
DO I = 1, 3
X(I,J) = -1.0 - j
ENDDO
ENDDO
DO I = 0, 2
H(I) = 1.0/(I+1)
ENDDO
CALL SFIRM1D(X(0,0), 1, 4, 0, 4, 5, H(0), 1, 0, 3,
& Y(0,0), 1, 6, 0, 6, ALPHA, BETA)
C/C++:
#include <scsl_fft.h>
float x[5][4], h[3], y[5][6];
float alpha = 1.0f;
float beat = 0.0f;
int i, j;
for (j=0; j<5; j++) {
Page 7
CFIRM1D(3S)CFIRM1D(3S)
x[j][0] = j + 1.0f;
for (i=1; i<4; i++) {
x[j][i] = -1.0f - j;
}
}
for (i=0; i<3; i++) {
h[i] = 1.0f/(i+1);
}
sfirm1d((float *) x, 1, 4, 0, 4, 5, h, 1, 0, 3,
(float *) y, 1, 6, 0, 6, alpha, beta);
The output is
Y(:,0) Y(:,1) Y(:,2) Y(:,3) Y(:,4)
Y(0,:) 1.0000 2.0000 3.0000 4.0000 5.0000
Y(1,:) -0.5000-1.0000-1.5000 -2.0000 -2.5000
Y(2,:) -1.1667-2.3333-3.5000 -4.6667 -5.8333
Y(3,:) -1.8333-3.6667-5.5000 -7.3333 -9.1667
Y(4,:) -0.8333-1.6667-2.5000 -3.3333 -4.1667
Y(5,:) -0.3333-0.6667-1.0000 -1.3333 -1.6667
Changing the values for ix0, ih0 and iy0 produces the following shifts in
the output:
ix0 = +1 Y(:,0) Y(:,1) Y(:,2) Y(:,3) Y(:,4)
Y(0,:) 0.0000 0.0000 0.0000 0.0000 0.0000
Y(1,:) 1.0000 2.0000 3.0000 4.0000 5.0000
Y(2,:) -0.5000-1.0000-1.5000 -2.0000 -2.5000
Y(3,:) -1.1667-2.3333-3.5000 -4.6667 -5.8333
Y(4,:) -1.8333-3.6667-5.5000 -7.3333 -9.1667
Y(5,:) -0.8333-1.6667-2.5000 -3.3333 -4.1667
ih0 = -1 Y(:,0) Y(:,1) Y(:,2) Y(:,3) Y(:,4)
Y(0,:) -0.5000-1.0000-1.5000 -2.0000 -2.5000
Y(1,:) -1.1667-2.3333-3.5000 -4.6667 -5.8333
Y(2,:) -1.8333-3.6667-5.5000 -7.3333 -9.1667
Y(3,:) -0.8333-1.6667-2.5000 -3.3333 -4.1667
Y(4,:) -0.3333-0.6667-1.0000 -1.3333 -1.6667
Y(5,:) 0.0000 0.0000 0.0000 0.0000 0.0000
iy0 = -1 Y(:,0) Y(:,1) Y(:,2) Y(:,3) Y(:,4)
Y(0,:) 0.0000 0.0000 0.0000 0.0000 0.0000
Y(1,:) 1.0000 2.0000 3.0000 4.0000 5.0000
Y(2,:) -0.5000-1.0000-1.5000 -2.0000 -2.5000
Y(3,:) -1.1667-2.3333-3.5000 -4.6667 -5.8333
Y(4,:) -1.8333-3.6667-5.5000 -7.3333 -9.1667
Y(5,:) -0.8333-1.6667-2.5000 -3.3333 -4.1667
Page 8
CFIRM1D(3S)CFIRM1D(3S)SEE ALSOCCOR1D(3S), INTRO_FFT(3S), INTRO_SCSL(3S
Page 9
