MODWT IMODWT code without fft

4 views (last 30 days)
Emiliano Rosso
Emiliano Rosso on 21 Mar 2023
Edited: Emiliano Rosso on 12 Apr 2023
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!

Answers (0)

Tags

Products


Release

R2020b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!