The filters used to produce the figures in the paper are available in FgrFltr.zip.
The available software consists of three independent packages.
1) Rational DWT and DFT-Modulated FBs
2) Filter Design
3) 2D Directional Transforms Derived From DFT-Modulated FBs
The zip file rational_dft.zip contains the MATLAB programs that implement the forward and inverse iterated orthonormal rational and overcomplete DFT-Modulated FBs along with a number of (almost) perfect reconstruction filter sets.
The iterated orthonormal rational FB consists of an iterated filter bank as shown below.
where p' = q-p. Let's load the filters for the case p = 5, q=6 and 4 vanishing moments.
>> load p5q6_vm4_set1;
Notice that the file names give the necessary structural information along with the number of vanishing moments.
>> who
Your variables are:
g_dft g_rat h_dft h_rat p q
'h_rat' and 'g_rat' are cell arrays holding the analysis and synthesis filter sets for the orthonormal rational FB.
>> h_rat
h_rat =
[1x65 double] [1x46 double]
h_rat{1} and h_rat{2} are the lowpass and highpass filters respectively. The filters in h_rat and g_rat are time reverses of each other.
>> max( abs( h_rat{1} - g_rat{1}(end:-1:1) ) )
ans =
0
Let's apply the iterated FB with 7 stages on a random signal of length 1000.
>> x = randn(1,1000);
>> no_stages = 7;
>> coef = DWT_pbyq(x,h_rat,p,q,no_stages)
coef =
[] [1x175 double]
[] [1x149 double]
[] [1x127 double]
[] [1x109 double]
[] [1x94 double]
[] [1x81 double]
[1x325 double] [1x71 double]
coef{n,k} holds the coefficients for the n^th stage, k^th band. Let's check perfect reconstruction.
>> y = IDWT_pbyq(coef,g_rat,p,q);
>> max( abs( x - y(1:length(x)) ) )
ans =
0.0010
We get almost perfect reconstruction. This can be verified by looking at the reconstruction error plotted along with the input signal via
>> plot(x);hold on; plot(x-y(1:length(x)),'r')
producing the figure below. The signal is plotted in blue and reconstruction error in red.
Now let's turn to overcomplete DFT-Modulated filter banks. Given a (p,q) pair, this is an FB as shown below (where W = exp(j2pi/q) ).
In the file we loaded, there are two more filters, 'h_dft' and 'g_dft'. These are the lowpass filters for the analysis and synthesis FB of the overcomplete DFT-Modulated FB. The remaining filters are obtained from these filters by frequency modulation. Again, these are time-reverse pairs.
>> max( abs( h_dft - g_dft(end:-1:1) ) )
ans =
0
Moreover, they are related to the lowpass filter of the rational FB by a normalization.
>> max( abs ( h_rat{1} - h_dft*sqrt(q) ) )
ans =
6.6613e-016
Let's decompose a random signal, this time using the iterated overcomplete rational DFT-Modulated DB with 9 stages.
>> x = randn(1,1000);
>> no_stages = 9;
>> c = DWT_DFT(x,h_dft,p,q,no_stages)
c =
[] [1x213 double] [1x213 double] [1x213 double] [1x213 double] [1x213 double]
[] [1x56 double] [1x56 double] [1x56 double] [1x56 double] [1x56 double]
[] [1x24 double] [1x24 double] [1x24 double] [1x24 double] [1x24 double]
[] [1x18 double] [1x18 double] [1x18 double] [1x18 double] [1x18 double]
[] [1x17 double] [1x17 double] [1x17 double] [1x17 double] [1x17 double]
[] [1x17 double] [1x17 double] [1x17 double] [1x17 double] [1x17 double]
[] [1x17 double] [1x17 double] [1x17 double] [1x17 double] [1x17 double]
[] [1x17 double] [1x17 double] [1x17 double] [1x17 double] [1x17 double]
[1x17 double] [1x17 double] [1x17 double] [1x17 double] [1x17 double] [1x17 double]
As in the transform above, c{n,k} holds the coefficients for the n^th stage, k^th band. Notice that there are q (=6) bands. Let's reconstruct the signal.
>> y = IDWT_DFT(c,g_dft,p,q);
>> max( abs( x-y(1:length(x) ) ) )
ans =
5.1006e-004
verifying the claimed almost perfect reconstruction. Since the transforms employ filters with complex coefficients, the reconstruction is complex as well. However, the imaginary part is negligible in this case as expected.
>> max( abs( imag( y ) ) )
ans =
1.0547e-014
The (almost) perfect reconstruction property can be verified by the following fragment of code as in the rational FB case
>> plot(x); hold on;
>> plot( x - real( y(1:length(x)) ), 'r' );
>> plot( imag( y(1:length(x)) ), 'y' );
which produces the figure below.
