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
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 toeplita 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, a can not find theoretical base for second step.
In the paper as references are 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:
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
|
DCTV |
DCTVII | ||||||
|
DSTVII |
DSTV | ||||||
|
DCTVI |
DCTVIII |
DCTII |
DCTIV | ||||
|
DSTVIII |
DSTVI |
DSTIV |
DSTII |
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: HOME Page , where you can read about these filters in more details.
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.
