Error during implementation of Wavelet based denoising using hybrid approach of adaptive flitering and single value decomposition

I want to denoise the signal using hybrid approach of adaptive flitering and single value decomposition. I have written the following code but it's showing error in Line 32

S_thresh = diag(S - threshold*(S > threshold));

Following is the code and signal image with nosie:

clear all

clc

load("Signal")

% Load the noisy signal or image

noisy_signal = C;

% Set parameters

wavelet = 'db2'; % Choose wavelet type

level = 5; % Decomposition level

threshold_type = 's'; % 's' for soft thresholding, 'h' for hard thresholding

threshold_multiplier = median(abs(C)) / 0.6745; % Threshold multiplier for adaptive thresholding

% Perform wavelet decomposition

[c, l] = wavedec2(noisy_signal, level, wavelet);

% Initialize denoised coefficients

denoised_c = cell(1, level + 1);

% Loop through each detail coefficient subband for denoising

for i = 1:level

% Extract detail coefficients for current subband

cd = detcoef2('all', c, l, i);

% Apply adaptive filtering

denoised_cd = zeros(size(cd));

for j = 1:size(cd, 1)

% Apply adaptive thresholding using singular value decomposition

[U, S, V] = svd(cd(j, :));

if threshold_type == 's'

threshold = threshold_multiplier * median(diag(S));

S_thresh = diag(S - threshold*(S > threshold));

elseif threshold_type == 'h'

S_thresh = diag(S.*(abs(S) > threshold_multiplier * median(diag(S))));

end

denoised_cd(j, :) = U * S_thresh * V';

end

% Update denoised coefficients

denoised_c{i} = denoised_cd;

end

% Keep approximation coefficients as is

denoised_c{level + 1} = appcoef2(c, l, wavelet, level);

% Reconstruct the denoised signal or image

denoised_signal = waverec2(denoised_c, l, wavelet);

% Display the results

figure;

subplot(1, 2, 1);

imshow(noisy_signal);

title('Noisy Signal');

subplot(1, 2, 2);

imshow(uint8(denoised_signal));

title('Denoised Signal');

[Pratyush Swain]

Hi Manindar,

It would be really helpful if you could share the exact error thrown during execution. Also sharing your signal data can help in reproducing the error you are facing.

Thanks

Answer 1

[MANINDER CHOUDHARY]

Hi pratyush,

Thank you for your feedback. Please find the link that contains signal that I intend to denoise using wavelet denoising that uses hybrid adaptive thresholding and single value decompositon (SVD) approach.Noise signal

The sampling frequeny is 1.5625e+08 Hz

Answer 2

[Pratyush Swain]

Hi Manindar,

The error in your script relating to "diag(S - threshold*(S > threshold))" is due to violation of the laws of Matrix Multiplication.

By carefully analyzing the variable sizes in your script, I observed both "threshold" and "S" have the same size of [1x2] hence "threshold x S" throws out an error. I suppose the value of "threshold" has to be a constant value,but it holds a matrix dimension due to use of "diag" function in line "threshold = threshold_multiplier * median(diag(S))". You can refer to dimension analogy as follows:

Variable Dimension

S -- [1X2]

diag(S) -- [2X2]

median(S) -- [1X2]

threshold -- [1X2]

Please note the behavior of diag function is as follows -

• diag(v) [here v is one dimensional vector]: returns a square diagonal matrix with the elements of vector v on the main diagonal.(i.e increases dimension)
• diag(A) [here A is a matrix]: returns a column vector of the main diagonal elements of A.(i.e decreases dimension)

Also using an inner "for-loop" for each row of detail coefficients matrix(cd) in line "for j = 1:size(cd, 1)" seems redundant as "svd" function can be applied on a matrix too.Please refer to workflow below:

for i = 1:level
    % Extract detail coefficients for current subband
    cd = detcoef2('all', c, l, i);

    % Apply adaptive thresholding using singular value decomposition
    [U, S, V] = svd(cd,"econ"); 
    % Used "econ" flag in the function above to work with large matrices
    
    % Perform thresholding as per threshold type - 'h'(hard) or 's'(soft)
    if threshold_type == 's'
        threshold = threshold_multiplier * median(diag(S));
        % Ommiting the use of extra diag function to maintain correct dimensions
        S_thresh = S - threshold*(S > threshold); 
    elseif threshold_type == 'h'
        S_thresh = S.*(abs(S) > threshold_multiplier * median(diag(S)));
    end
    
    denoised_cd = U * S_thresh * V';
    % Update denoised coefficients
    denoised_c{i} = denoised_cd;
    
end

Please refer to following documentations for more information:

1. https://www.mathworks.com/help/matlab/ref/double.svd.html

2. https://www.mathworks.com/help/matlab/ref/diag.html

3. https://www.mathworks.com/help/matlab/ref/mtimes.html

Hope this helps you in progressing ahead.

[MANINDER CHOUDHARY]

Dear Pratyush,

Thank you for the detailed explaination and your valuable input. The error still continues. The reconsuction line shown an error. I'm not able to figure out why. I also tired to remvoe 2 as my signal is 1-D so i used commmand waverec rather waverec2 or appcoeff2

% Reconstruct the denoised signal,

denoised_signal = waverec(denoised_c, l, wavelet);

[Pratyush Swain]

Hi Maninder,

The issue of reconstruction is another task in itself. The waverec function expects "denoised_c" as a single vector of wavelet coefficients but currently its in form of an cell array containing only thresholded detail coefficients at varied levels.

So a method/logic has to be figured out how to reconstruct the wavelet coefficients. This thread might prove helpful: https://www.mathworks.com/matlabcentral/answers/213354-how-can-i-combine-detail-coefficients-from-different-levels-of-a-wavelet-decomposition

Thanks

[MANINDER CHOUDHARY]

I'm trying to use denoise the signal using Dual tree complex variable transform but there is no library in Matlab. Can you please guide how to denoise this signal using Dual treee complex wavelet transform?