Implementation of Fourier Domain Denoising with Hard Threshold

I just tried the Fourier denoising method with a hard threshold and my code is as follow:

F = imnoise(phantom(128),'gaussian');
B = fft(F);
s = 0.3;
for i = 1:128
for j = 1:128
if abs(B(i,j))<s
B(i,j)=0;
end
end
end

C=ifft(B);
subplot(121)
imshow(F,[])
subplot(122)
imshow(C,[])


Left is the phantom after adding noise and the right is the denoised image from my code.

The Fourier denoising hard threshold method just uses threshold value to keep high frequency coefficients and the coefficients below the threshold to be 0.

In my code, I just added gaussian noise with default mean 0 and variance 0.01 to my phantom image and then I set the threshold to be 3 times of the standard deviation of the noise, which is 0.3 (this is some good practical threshold I looked at from literature). But why the result is not good? It seems it is still very noisy. Even I change my threshold the result is not good. Why? And I want to my denoised image as close to the original phantom.

• You should use fft2 (2-D) instead of fft (1-D). Jun 9, 2015 at 10:37
• Have you checked the magnitudes of your FFT? I don't think fft or fft2 return normalized results. 0.3 might be good for the spatial domain noise, but be far too small for the transform. Someone correct me if I'm wrong. By the way, the resulting image is identical to the input except for a vertical band trimmed from the middle. Nov 3, 2015 at 19:30
• SE.DSP wishes you a happy new year 2017, with a kind reminding signal that your question or its answers may require some action from you (edit, update, votes, acceptance, etc.) Jan 15, 2017 at 16:24

Thresholding in the Fourier domain is an archaic method sometimes called spectral subtraction, used for background noise removal in speech. Bad results in image processing can be due to several factors:

1. Misinterpretation: you say "The Fourier denoising hard threshold method just uses threshold value to keep high frequency coefficients". Not quite. You keep "high amplitude coefficients", no matter low or high frequency.
2. Test data: the Shepp-Logan model is known to be quite specific. Bad results on it do not imply bad results on real data, yet:
3. Transform choice: Fourier is very good for close-to-stationnary data. Standard images with smooth parts inside contours are not well represented in the Fourier domain, at least not concentrated enough for acurate data compaction and reconstruction with few (high) coeffcients. Derivatives, or multiscale directional derivatives like wavelets can help better.
4. Implementation: Fourier implementation via FFT suffers boundary artifacts coefficients, not talking about the 1D implementation you seem to use (prefer fft2 here)
5. Coefficient selection: hard and scalar thresholding typically yields more important visual artifacts than "smoother" thresholding and vector or group thresholding.
6. Formulation: you can generally get better results when you do not use a "parameter-tuned algorithm", but prefer a variational formulation based on your priors, including noise statistics.

So you have plenty of room to improve your results. Two references for multivariate image denoising with automated threshold (for Gaussian noise):

You added noise with Variance of 0.01 in the spatial domain.
Now you need to calculate what's the variance of the Noise per each bin in the DFT.

It's not the same, there are gains in the DFT domain, otherwise no one would use it.

By the way, there are much faster methods to apply "Hard Threshold" than the loop in your program. Use logical indexing to do it faster.

• But the variance in the fourier domain, as a matter of fact, is just a different number. I have tried all kinds of numbers as a threshold, it just did not work well...
– Ian
Jun 9, 2015 at 15:45
• What do you mean when you say "It just did not work well..."? Do you mean you lost edges and other data?
– Royi
Jun 9, 2015 at 17:05
• @Ian, Could you address what's missing in my answer in order to mark it?
– Royi
Jun 30, 2022 at 8:43

You should use the 2D Fourier transform function fft2 instead of fft. When provided with a 2D array arguement, fft returns 1D FFT results for the columns.

In order to remove the noise you should apply a low pass filter, as opposed to zeroing coefficients with small magnitudes. By using the fftshift and meshgrid functions you can easily do this in frequency space.

• Thanks. But why can't I zero coefficients with small magnitudes because the noise are supposed to have low magnitudes and why the denoising result is not good, as the images I showed?
– Ian
Jun 10, 2015 at 2:49

In literature it is shown that direct thresholding does not perform well. Consider iterative thresholding. You can find some papers on that ( see basis pursuit, matching pursuit etc.,). You can also read this paper - http://arxiv.org/abs/1504.00976