ti_denoise.m

用maltlab编程实现了平移不变量去噪的快速算法

function y=TI_Denoise(signal,wavename,L)

%  Fast TI_Denoising of 1-d signal with wavelet thresholding.
%  Usage 
%    y=TI_Denoise(signal,wavename,L)
%  Inputs
%    signal	    1-d Noisy signal, length(Noisy)= 2^J. J must be an positive integer.
%    wavename   name of wavelet   
%    L      	Low-Frequency cutoff for shrinkage. L << J
%  Outputs 
%    y     	    the signal after being denoised



%%%Fast TI forward wavelet transform

n=length(signal);
%Initialize TI table
table=zeros(n,L+1);
table(:,1)=signal';
%Different decompose scales from 1 to L
for i=1:L,
    %Size of the low_frequency coefficients or wavelet coefficients at different scales
    size=n/2^i;
    %the low_frequency coefficients of the prior scale
    low_coef=table(:,1);
    for j=0:2:2^i-2,
        %Get the low_frequency coefficients belta(J+1-i,j/2) 
        sig0=low_coef(size*j+1:size*(j+2));
        %DWT to get the low_frequency coefficients belta(J-i,j) and wavelet coefficients alpha(J-i,j)
        [a1,b1]=dwt(sig0,wavename,'mode','per');
        %Store the coefficients into TI table
        table((j*size+1):(j+1)*size,1)=a1;
        table((j*size+1):(j+1)*size,i+1)=b1;
        %Shift 1 on belta(J+1-i,j/2) 
        sig1=shift_left(sig0);
        %DWT to get the low_frequency coefficients belta(J-i,j+1) and wavelet coefficients alpha(J-i,j+1)
        [a2,b2]=dwt(sig1,wavename,'mode','per');
        %Store the coefficients into TI table
        table(((j+1)*size+1):(j+2)*size,1)=a2;
        table(((j+1)*size+1):(j+2)*size,i+1)=b2;
    end
end


%%%Soft threshold denoising

%estimate the standard deviation of additive noise.
d=table(:,2);
sigma=median(abs(d))/0.6745;
%Set the threshold
thresh=sqrt(2*log(n))*sigma;
%Soft threshold
for i=2:L+1,
    for j=1:n,
        if abs(table(j,i))<thresh,
            table(j,i)=0;
        elseif table(j,i)>=thresh,
            table(j,i)=table(j,i)-thresh;
        else
            table(j,i)=table(j,i)+thresh;
        end
    end
end


%%%Reconstruct the signal

%Different reconstruction scales from L to 1
for i=L:-1:1,
    size=n/2^i;
    for j=0:2:2^i-2,
        %get the low_frequency coefficients belta(J-i,j) and wavelet coefficients alpha(J-i,j)
        a1=table((j*size+1):(j+1)*size,1);
        b1=table((j*size+1):(j+1)*size,i+1);
        %IDWT 
        sig0=idwt(a1,b1,wavename,'mode','per');
        %get the low_frequency coefficients belta(J-i,j+1) and wavelet coefficients alpha(J-i,j+1)
        a2=table(((j+1)*size+1):(j+2)*size,1);
        b2=table(((j+1)*size+1):(j+2)*size,i+1);
        %IDWT
        sig1=idwt(a2,b2,wavename,'mode','per');
        %Inversely shift 1
        sig1=shift_right(sig1);
        %Get the low_frequency coefficients belta(J+1-i,j/2) 
        low_coef(size*j+1:size*(j+2))=(sig0+sig1)/2;
    end
    %Update the low_frequency coefficients
    table(:,1)=low_coef;
end
%return the denoised signal
y=table(:,1);