Fortran code for signal processing

Signal processing is a domain where Fortran was once dominant but has now been largely replaced by C, C++, and Matlab. R has a package gsignal with many contributors. It calls C++ code that looks like Fortran 77 with braces. The NumericMatrix class used accesses elements as x(i,j). I wonder how efficient it is compared to the native arrays of Fortran. If it could be shown that for some algorithms, Fortran code is both shorter and faster than a C++ implementation, maybe some of the computing done in R would be done in Fortran, and maybe more R packages would be written in modern Fortran. The C++ code in the package is not the code one would use when performance is critical, but ideally one could show that performance-critical Fortran code is fast, easy to read, and safe.

3 Likes

In my university library I have seen several older textbooks on signal processing with Fortran.

Digital filter design

Author: Parks, Thomas W.
Place of Publication, Publisher, Year of Publication: New York u.a., Wiley, 1987

Appendix contains Fortran program for FIR and IIR filters.

Adaptive digital filters and signal analysis

Author: Bellanger, Maurice
Place of Publication, Publisher, Year of Publication: New York u.a., Dekker, 1987

Some Fortran procedures are mentioned in the table of contents.

Digital signal processing in Fortran

Author: Taylor, Fred
Place of Publication, Publisher, Year of Publication: Lexington, Mass. [u.a.], Lexington Books, 1976

Digital signal analysis

Author: Stearns, Samuel D.
Edition: 2. ed.
Place of Publication, Publisher, Year of Publication: Englewood Cliffs, NJ, Prentice Hall, 1990

Appendix B provides algorithms in Fortran-77

Practical Applications in Digital Signal Processing

Author: Newbold, Richard
Edition: 1st edition
Place of Publication, Publisher, Year of Publication: [Erscheinungsort nicht ermittelbar], Pearson, 2012

Appendix A is titles “Mixed Language C/C++ FORTRAN Programming”. Even though the book is published in 2012, it appears the authors were unaware of the standardized Fortran/C interoperability and rely on some manual mangling/un-mangling protocol. Section A.10 has the following nice paragraph:

The Parks-McClellan FORTRAN program for computing optimum linear phase digital filters has been around in one form of legacy code or another since December 1973. It can be found in practically every signal processing engineer’s toolbox. It’s probably in the possession of a huge majority of engineers who were actively involved in designing digital filters back in the 1970s and 1980s. Most of the glossy, Windows-based filter design programs available today use either the Parks-McClellan program or some variant of it.

Quite a testament! The original Parks-McClellan program can be found in:

McClellan, J., Parks, T. W., & Rabiner, L. (1973). A computer program for designing optimum FIR linear phase digital filters. IEEE Transactions on Audio and Electroacoustics , 21 (6), 506-526. https://doi.org/10.1109/TAU.1973.1162525

The main author is an associate professor in behavioral sciences and neuropsychology. I suspect the many contributors might be coming from how CRAN assigns contributions if sub-packages are called. The R document Who Did What? states contributors do not appear in the package citation.

The gsignal R package is a reimplementation of the octave Signal processing package. I’ve had a look at the source code, and there is no Fortran to be found there. Only the underlying algorithms in C++ with the MATLAB wrapper routines.

The signal processing I know should include time series analysis and FFT.
We may be able to build related signal processing packages based on rafat/ctsa and fftpack, or something.

1 Like

NumericMatrix is from the package Rcpp. I have used it extensively. Passing a large matrix from R to C++ using NumericMatrix can avoid copying, which makes it faster than passing the matrix from R to Fortran using R’s built-in Fortran FFI. Luckily, the R package dotCall64 has solved this problem. My personal experience is that Fortran is as fast as C++ when dotCall64 is used.