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 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:
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:
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.
| Illustration demo overlap-add method |
