Digital Signal processing algorithms

In this page I will place description and implementation of various DSP algorithms. At least Matlab realization will be presented. Implementation of some algorithms can be presented also on various programming languages. This page is now without some ordering system. The necessity of ordering comes with increasing of number of algorithms.

Solving toeplitz system y=Px
R. Kumar. A fast algorithm for solving a toeplitz system of equations. IEEE Trans. on Acous., Speech, and Sign. Proc., 33:254--267, February, 1985

MatLab implementation

This algorithm solves toeplitz equation y=Px, where P is toeplitz matrix. Such kind of equations arise for example in inverse filtering problems Computational complexity of this algorithm O(n*log2(n)^2) is less than complexity of well-known Lewinson-Durbin algorithm O(n^2).

Algorithm consists of three steps
Step 1:
- Constructing circulant matrix based on original toeplitz matrix and inverting obtained circulant matrix by means of FFT.
Step 2:
- First row of inverse matrix P^(-1) is found by means of polynomial division algorithm.
Step 3:
- Vector 'x' is found. Here is used the first row of inverse matrix and special structure of original matrix P.

Although this algorithm works, I can not find theoretical base for second step.This step references to next papers
 W.F. Trench. "An algorithm for the inversion of finite Toeplitz matrices", J.SIAM, vol.12, pp.512-522, Sept. 1964.
 S. Zohar. "The solution of a Toeplitz set of linear equation", J. Ass. Comput. Mach., vol.21, pp.272-276, Apr. 1974.
 H.H. Rosenbrock. State Space and Multivariable Theory. New-York: Wiley, 1970.

FIR filtering using filter banks
M.Vetterli. Running FIR and IIR Filtering Using Multirate techniques. IEEE Trans. On Acoustics, Speech and Signal Processing, 36:730-738, May 1988.

Good MatLab implementation you can find in Jelen Tesic's report

Short description:
It is well-known that by filter length less than 100, spectral methods for computing convolutions does not have any benefits in comparison with direct methods. In this case the best way is to use algorithms based on chinese remainder theorem . However, for length of several tens of samples, obtained algorithms (nested algorithms) are sophisticated. Additionally by changing of filter length you need to synthesize another algorithm - for new filter length.
Idea of presented algorithm is based on separating impulse response of the filter and input signal on polyphase components. In this case convolution can be presented as pseudocyclic convolution of polyphase components of the filter and input signal. Fast algorithm is synthesized by means of chinese remainder theorem.
Benefits of the algorithm:
1. Convolution is divided on several shorter convolutions. Each convolution is computed on reduced sampling frequency.
2. Whole amount of computation can be reduced on 25 and more percents.
Drawbacks:
1. Bigger code in comparison to code realizing direct calculation.
2. Bigger data size.

Afterword:
1. In original paper of Martin Vetterli there is no any reference to chinese remainder theorem, so I made it in papers Eigen Transforms over Finite Rings in Fiter Bank Structures and New Approach to Realizing Digital Filters Using Filter Banks and Eigen Transforms in Polynomial Rings. In these paper you can also find extension of this approach on the case of finite fields.
2. Only at spring 2005 I found occasionally, that in two years before me, Prof. Mitra in his paper Ing-Song Lin, and Sanjit K. Mitra, "Overlapped Block Digital Filtering", IEEE TRANSACTIONS ON CIRCUITS AND SYSTEMS-II, VOL. 43, NO. 8, AUGUST 1996 discovered this fact too. So I should accept his priority in discovering these algorithms. I can only say, that at the time when I wrote my papers, I had no access to latest IEEE journals.

Fast algorithm of 50-point Fourier transform

“In fact good FFT algorithms exist for any block length”. (Richard Blahut)

Fast realization of discrete Fourier transform of power of two length by means of Cooley-Tukey approach is well-known. Of course this algorithm is most effective and most uniform structured. (Often I choose transforms of power of two length too). However, 50-point Fourier transform algorithm illustrates above depicted Blahut’s expression very well. By proper using of several approaches described in Blahut's book "Fast Algorithms for Digital Signal Processing" for splitting large-point Fourier transforms into smaller ones computational complexity of algorithm for every length can be significantly decreased.

Combinatorial square root algorithm
 Algorithm description Graph of algorithm

I never meant that square root can be computed so quick and so easy. Result is computed by means of n/2 iteration, where n is number of bits for representing original number. In each iteration only one subtraction and one comparison operations are used. The rest are shifts and concatenation operations. That's why this algorithm is well suited for FPGA realization.

Similar algorithms are developed for many other mathematical functions.
I recommend good books on this topic:

• Computer Arithmetic Algorithms 2nd edition, by Israel Koren, Natick, MA: A. K. Peters, 2002
• Muller, J. M., Elementary functions: algorithms and implementation, Birkhauser, 1997.

I would like to point out very interesting paper
Auger, F. Zhen Lou Feuvrie, B. Feng Li. Multiplier-Free Divide, Square Root, and Log Algorithms [DSP Tips and Tricks]
There are matlab files illustrating all algorithms presented in this paper.

Discrete Trigonometric Transform
Fast methods for cosine and sine transforms computation through one another
 Markus Puschel and Jose Moura The Algebraic Approach to the Discrete Cosine and Sine Transforms and their Fast Algorithms SIAM Journal of Computing 2003, Vol. 32, No. 5, pp. 1280-1316 R. Gluth, “ Regular FFT-related transform kernels for DCT/DST- based polyphase filter banks,” Proc. IEEE ICASSP 1991, pp.2205-2208, Toronto, Canada, May 1991. Mathematical relations between discrete trigonometric transform (file in Word-format) Matlab code presenting some transform relations

There are 8 types of discrete cosine and sine transforms. In a paper the first link points Markus Puschel and Jose Moura were developed beautiful theory that combines all 16 trigonometric transforms in one framework and shows that sine and cosine transforms of types 1...4 are related to one another in a very simple fashion. In the same fashion are related to one another sine and cosine transforms of types 5...8. Second link consist paper of Rolf Gluth in that effective approach for computation of sine and cosine transform is presented. Hence, sine and cosine transforms of types 1...4 and Fourier transform (9 most useful transforms) are related to one another in strong mathematical relations. It means that if we have in our DSP-library the code for only one transform (offen only FFT code is available), we can compute any other of 8 transforms using existing code by the small cost of additional pre- and postcalculation. Exhaustive information about this topic you can get visiting homepage of Markus Puschel.

Relations between Discrete Trigonometric Transforms
Matlab reports

• May be it is interesting for someone - I placed here matlab scripts GenDTT.zip for generating transform generations. Just place proper pair of transforms at the line 143 and run. Just be sure that transform pair is from cells of the same color.
• I made detailed comments to the DST2<->DCT1 transform conversion example presented in the M. Pueschel's paper mentioned above (DTT transform relationship derivation(detailed).pdf).

EMQF filters
 Implementation demo for EMQF filters Matlab files used in designing EMQF filters

EMQF filters (elliptic filters with minimal Q-factor) are particular case of elliptic filters that possess the same ripples in stop- and passband. Despite of this restriction, this class of filters have very interesting properties on the base of these we can design effective IIR filters, based on coupled allpass structure, we can realize effective IIR decimators and interpolators, Hilbert pair IIR filters etc. Beside demo page for EMQF you can visit Miroslav D. Lutovac: publication list , where you can read about these filters in more details.

# Design Examples Using Least Square IIR Filter Optimization Program

 Implementation demo for IIR filter design Matlab files used in designing IIR filters I recreated all matlab scripts from Mr. Lang's dissertation as well. You can find it under the link diss_lang.zip.

This package is used for IIR filter design using prescribed magnitude and phase responses using least square optimization. One of the parameters is maximal radii of filter poles. It allows guaranteeing filter stability by coefficient discretization and controlling dynamic range of internal filter state. Package is borrowed from dissertation Vienna, June 1999.

# Conjugate Symmetric Sequency-Ordered Complex Hadamard Transform

 Illustration demo on CS-SCHT Matlab package for using CS-SCHT

This package was written after the paper "Conjugate Symmetric Sequency-Ordered Complex Hadamard Transform", Aye Aung   Boon Poh Ng   Rahardja, S.. Authors claimed that fast algorithm for real type of this transform takes about 25% less operations than conventional Hadamard transform. Basis functions of this transform are approximations of DFT basis functions that makes this transform attractive tool for quick (though inaccurate) spectrum estimation. In illustration demo you can see application of this transform to image compression and compare performance of it using different transforms.

# Filter visualization tool (fvtool) for Octave

 Illustration demo on fvtool for Octave Matlab script (fvtool.m and demo-script)

Sometimes on my mobile devices I'm working in Octave (open source equivalent of Matlab). Absence of filter visualization tool (fvtool) was really inconvenient. So at one day I created fvtool script for Octave. It is not so pretty and powerful as Matlab tool (in Matlab it is interactive), but provide major functionalities of original fvtool from Matlab. You can draw various responses: frequency, phase, phase delay, group delay responses and impulse response as well. Frequency and magnitude axes can be displayed in a linear or logarithmic scale. For convenience of analysis, samplerate can be set i order to provide correct displaying of axis values.

# Overlap-Add and Overlap-Save Block Convolution Methods Using FFT of Size N

 Illustration demo overlap-add method Illustration demo overlap-save method Matlab scripts for overlap-add and overlap-save methods

Classic overlap-add and overlap-save methods are using FFT size of 2N, where N is block size. Proposed methods provide reduction of calculation amount by replacing FFT size of 2N by so named skew transform of size N that is an eigentransform related to convolution mod(x^N-j). The above mentioned transform is a skew transform and can be easily implemented as multiplication by diagonal matrix of twiddle coefficients exp(j*pi/(2N)*k), k=0..N-1, and appplying FFT of size N. Regarding this consideration computational amount is reduced significantly. Original idea was proposed in the paper Jung Gap Kuk, Seyun Kim, and Nam Ik Cho, "A New Overlap Save Algorithm for Fast Block Convolution and Its Implementation Using FFT," Journal of Signal Processing Systems(Online), Feb 2010.