MODWT IMODWT code without fft
5 views (last 30 days)
Show older comments
Hi.
I studied the maximal overlap wavelet transform and its properties on "Wavelet Methods for Time Series Analysis by Donald B. Percival, Andrew T. Walden " and I saw that the application of the fft is performed only after the development of the algorithm for speed up the code. To better understand how it works I need to develop this algorithm without using the fft. Is there any library that does this?
Thank you!
Now i found the modwt code :
function [coeffs] = NOFFTmodwt(x, filters, level)
% filters: 'db1', 'db2', 'haar', ...
% return: see matlab
% filter
N=length(x);
[Lo_D,Hi_D] = wfilters(filters);
g_t = Lo_D / sqrt(2);
h_t = Hi_D / sqrt(2);
wavecoeff= cell(1,5);
v_j_1 = x;
for j = 1:level
w = circular_convolve_d(h_t, v_j_1, j);
v_j_1 = circular_convolve_d(g_t, v_j_1, j);
wavecoeff{j} = w;
end
wavecoeff{level+1} = v_j_1;
coeffs = vertcat(wavecoeff{:});
coeffs(1:level,1:N)=-coeffs(1:level,1:N);
function w_j = circular_convolve_d(h_t, v_j_1, j)
% jth level decomposition
% h_t: \tilde{h} = h / sqrt(2)
% v_j_1: v_{j-1}, the (j-1)th scale coefficients
% return: w_j (or v_j)
N = length(v_j_1);
L = length(h_t);
w_j = zeros(1,N);
l = 0:L-1;
for t = 1:N
index = mod(t -1 - 2^(j-1)*l, N) + 1;
v_p = v_j_1(index);
w_j(t) = dot(h_t, v_p);
end
It's freely translated from the Python function modwtpy which you can find here :
It's equal in results to the Matlab modwt function , border sx oriented for the 'haar' , shit invariante , same coeffs, ecc..
Now I must find the same for the imodwt. this is what i managed to do but it doesn't work!
function [x] = NOFFTimodwt(coeffs, filters)
% inverse modwt
% filter
[Lo_R,Hi_R] = wfilters(filters);
g_t = Lo_R / sqrt(2);
h_t = Hi_R / sqrt(2);
level = size(coeffs, 1) - 1;
v_j = coeffs(end, :);
for j = level:-1:1
%w_j = circshift(coeffs(j,:),[0,2^(j-1)]);
%v_j = circular_convolve_s(h_t, g_t, w_j, v_j, j);
v_j = circular_convolve_s(h_t, g_t, coeffs(j,:), v_j, j);
end
x = v_j;
function v_j_1 = circular_convolve_s(h_t, g_t, w_j, v_j, j)
% (j-1)th level synthesis from w_j, w_j
% see function circular_convolve_d
N = length(v_j);
L = length(h_t);
v_j_1 = zeros(1,N);
l = 0:L-1;
for t = 1:N
index = mod(t -1 + 2^(j-1)*l, N) + 1;
%index = mod(t + 2^(j-1)*l - 1, N) + 1;
w_p = w_j(index);
v_p = v_j(index);
v_j_1(t) = dot(h_t, w_p) + dot(g_t, v_p);
end
Any suggestions?
Thanks!
0 Comments
Answers (0)
See Also
Categories
Find more on Discrete Multiresolution Analysis in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!