MODWT IMODWT code without fft

2 ビュー (過去 30 日間)
Emiliano Rosso
Emiliano Rosso 2023 年 3 月 21 日
編集済み: Emiliano Rosso 2023 年 4 月 12 日
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 件)

カテゴリ

Help Center および File ExchangeDiscrete Multiresolution Analysis についてさらに検索

タグ

製品


リリース

R2020b

Community Treasure Hunt

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

Start Hunting!

Translated by