Numerical differentiation based on wavelet transforms

Jianwen Luo (view profile)

14 Feb 2007 (Updated )

Numerical differentiation based on wavelet transforms (CWT and DWT)

dudx=derivative_dwt(u,wt_name,wt_level,dx,trt_flag)
```function dudx=derivative_dwt(u,wt_name,wt_level,dx,trt_flag)
%
% Differentiation (Derivative) of Sampled Data Based on Discrete Wavelet Transform
%
% dudx=derivative_dwt(u,wt_name,wt_level,dx)
%
% u:        uniformly-sampled data
% wt_name:  name of the wavelet function (haar or spl)
% wt_level: level of wavelet decomposition
% dx:       sampling interval (default=1)
% trt_flag: flag of translation-rotation transformation for boundary effect (default=1)
% dudx:     differentiations of data (u)
%
%           derivative_cwt
%
% Reference:
% J. W. Luo, J. Bai, and J. H. Shao,
% "Application of the wavelet transforms on axial strain calculation in ultrasound elastography,"
% Prog. Nat. Sci., vol. 16, no. 9, pp. 942-947, 2006.

if nargin<5
trt_flag=1;
end
if nargin<4
dx=1;
end

if trt_flag
x=(1:length(u))*dx;
a=(u(end)-u(1))/(x(end)-x(1));
b=u(1)-a*x(1);
u=u-a*x-b;
else
a=0;
end

wt_name=lower(wt_name);

if strcmp(wt_name,'haar')
h0=[sqrt(2)/2 sqrt(2)/2];  %the decomposition low-pass filter
h1=[-sqrt(2)/2 sqrt(2)/2]; %the decomposition high-pass filter
elseif strcmp(wt_name,'spl')
h0=[0.125 0.375 0.375 0.125]*sqrt(2);
h1=[-2 2]*sqrt(2);
else
error('wavelet name error');
end

y0=u;

% Algorithme a Trous
for n=1:wt_level
h0_atrous=[h0' zeros(length(h0),2^(n-1)-1)]';
h0_atrous=h0_atrous(1:(length(h0)-1)*(2^(n-1)-1)+length(h0));

h1_atrous=[h1' zeros(length(h1),2^(n-1)-1)]';
h1_atrous=h1_atrous(1:(length(h1)-1)*(2^(n-1)-1)+length(h1));

y1=conv(y0,h1_atrous);
y0=conv(y0,h0_atrous);
end

index=round(length(y1)/2-length(u)/2)+[1:length(u)];
dudx=y1(index);

wt_scale=2^wt_level;

if strcmp(wt_name,'haar')
dudx=-dudx/wt_scale^(3/2)*4;
elseif strcmp(wt_name,'spl')
dudx=-dudx/wt_scale^(3/2);
else
error('wavelet name error');
end

dudx=dudx/dx+a;

```