Matlab – Low pass gaussian filter with a specified cut off frequency

gaussianimage processingMATLAB

I'm playing around with hybrid images, and wanted to use a gaussian filter to low pass filter an image. However, to make hybrid images, 2 filters are supposed to be used on the 2 images being combined with different cut off frequencies.

Does fspecial() allow us to specify cut off frequencies when we employ it to make a gaussian filter? (I know that we can specify the filter size and sigma and that there is some relation between sigma and cut off frequency). If we can only specify cut off frequencies using sigma, then what sigma would I need to set to get a cut off frequency of say 0.2 Hz.

Best Answer

I'll first answer regarding 1D and the rest will follow. It may look trivial, but bear with me for a while. Lets assume the following code:

t=linspace(0,20,2^10); %time vector in seconds
w=0.2; %in Hz
signal=cos(2*pi*w*t)+rand(1,length(t))-0.5; % signal in seconds
dt=t(2)-t(1) ;
N=length(signal);
df=1/(N*dt);    % the frequency resolution (df=1/max_T)
if mod(N,2)==0
    f_vec= df*((1:N)-1-N/2);  % for EVEN length vectors
else
    f_vec= df*((1:N)-0.5-N/2); 
end

So, we have created a noisy signal of a specific frequency. f_vec is the frequency vector that stretches from f =[-f_max,-f_max+df,...,0,...,f_max], with f_max=1/(2*dt). If we now design a 1D Gaussian filter (in the fourier space) as follows: 

f_vec0=0;
sigma=1;
filter=exp( -(f_vec-f_vec0).^2./(2*sigma^2));

and then filtering in the fourier doamin:

f_signal=fftshift(fft(signal));
filt_signal=fftshift(ifft(f_signal.*filter));

So, from the filter we applied, sigma=1 means the the cut off frequency (that I decided is 1% of the filter's maximum (which is 1)) is approximately at 3 Hz:

 cutoff_freq=f_vec(find(filter>=0.01,1,'last'))

Taking this to 2D is trivial, just be careful with units. For images, there are pixels as the position unit, and 1/pixels as the spacial frequency. The fspecial function generates a 2D matrix of a predefined filter. The usage of fspecial is usually like this:

PSF = fspecial('gaussian',hsize,sigma);
Blurred = imfilter(Image,PSF,'symmetric','conv');

Using convolution is just like multiplying in the Fourier domain. The sigma in the Fourier domain is proportional to 1/sigma of the position domain, etc...

Related Topic