Vadim Kudryavtsev home page

** **

**
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
- FIR filtering using filter banks
- Fast algorithm of 50-point Fourier transform
- Combinatorial square root algorithm
- Discrete Trigonometric Transform
- EMQF filters
- Design Examples Using Least Square IIR Filter Optimization Program
- Conjugate Symmetric Sequency-Ordered Complex Hadamard Transform
- Filter visualization tool (fvtool) for Octave
- Overlap-Add and Overlap-Save Block Convolution Methods Using FFT of Size N

**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

Page 1-2 Page 3-4 Page 5-6 Page 7-8 Page 9-10 Page 11-12 Page 13-14

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**.

**Comments:**

Although this algorithm works, I can not find theoretical base for
second step.This step references to next papers

[4] W.F. Trench. "An algorithm for the inversion of finite Toeplitz
matrices", J.SIAM, vol.12, pp.512-522, Sept. 1964.

[5] S. Zohar. "The solution of a Toeplitz set of linear equation", J.
Ass. Comput. Mach., vol.21, pp.272-276, Apr. 1974.

[15] 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.

Page 1-2 Page 3-4 Page 5-6 Page 7-8 Page 9-10

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 [1].
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**

Description of algorithm developing

Algorithm realization on Matlab

“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

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

Additional notes:

- 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.

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 Mathias Lang, "Algorithms for the Constrained Design of Digital Filters with Arbitrary Magnitude and Phase Responses", Vienna, June 1999.

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.

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.

Illustration demo overlap-add method |

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.