The zip file rational_oc.zip contains the MATLAB programs that implement the forward and inverse overcomplete rational DWT and some perfect reconstruction filter sets for different rational sampling factors.
The overcomplete rational DWT consists of an iterated filter bank as shown in the figure below.
Let us demonstrate how to use the code using a filter set for the p=3, q=4 case. First let's load the filters.
>> load Q4_K4_N5
The file name indicates the properties of the filters in the filter set. The number after Q indicates the 'q' in the above figure. For the included filter sets 'p' is taken to be equal to 'q-1'. The number after K indicates the number of vanishing moments for the highpass filters. The number after N indicates the number of regularity factors of the lowpass filter.
>> who
Your variables are:
g h
'h' and 'g' are cell arrays holding the filters for the forward and inverse transforms. For this setting (p=3,q=4), we need 5 (i.e. q+1) filters.
>>
h =
[1x29 double] [1x10 double] [1x10 double] [1x11 double] [1x13 double]
Since this is a tight frame, the inverse filters are time reverses of the forward filters.
>> max( abs( h{2} - g{2}(end:-1:1) ) )
ans =
0
We will denostrate the use of the forward transform with 8 stages. We apply the transform on a random signal of length 1000.
>> p = 3; q = 4; no_stage = 8;
>> x = randn(1,1000);
>> coef = DWT_pbyq_oc(x,h,p,q,no_stage)
coef =
[] [1x253 double] [1x253 double] [1x253 double] [1x253 double]
[] [1x192 double] [1x192 double] [1x192 double] [1x193 double]
[] [1x146 double] [1x146 double] [1x147 double] [1x147 double]
[] [1x112 double] [1x112 double] [1x112 double] [1x113 double]
[] [1x86 double] [1x86 double] [1x87 double] [1x87 double]
[] [1x67 double] [1x67 double] [1x67 double] [1x68 double]
[] [1x53 double] [1x53 double] [1x53 double] [1x53 double]
[1x125 double] [1x42 double] [1x42 double] [1x42 double] [1x43 double]
coef{n,k} holds the coefficients for the n^th stage, k^th band. Let's check perfect reconstruction.
>>y = IDWT_pbyq_oc(coef,g,p,q);
>> y = y(1:length(x));
>> max(abs(x-y))
ans =
1.1657e-15
We can also look at the discrete time atoms at any stage using the inverse transform.
>> x = zeros(1,50);
>> coef = DWT_pbyq_oc(x,h,p,q,no_stage);
>> coef{8,2}(end) = 1;
>> y = IDWT_pbyq_oc(coef,g,p,q);
>> plot(y)
This fragment of code produces the figure below.
