Clear Filters
Clear Filters

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

4 views (last 30 days)
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');
  1 Comment
Pratyush Swain
Pratyush Swain on 15 May 2024
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

Sign in to comment.

Answers (2)

MANINDER CHOUDHARY
MANINDER CHOUDHARY on 15 May 2024
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

Pratyush Swain
Pratyush Swain on 15 May 2024
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:
Hope this helps you in progressing ahead.
  3 Comments
Pratyush Swain
Pratyush Swain on 16 May 2024
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
MANINDER CHOUDHARY on 23 May 2024 at 13:46
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?

Sign in to comment.

Categories

Find more on Denoising and Compression in Help Center and File Exchange

Products


Release

R2024a

Community Treasure Hunt

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

Start Hunting!