我将首先回答有关一维的问题,其余的将随之而来。它可能看起来微不足道,但请耐心等待一段时间。让我们假设以下代码:
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
因此,我们创建了特定频率的噪声信号。 f_vec 是从 f =[-f_max,-f_max+df,...,0,...,f_max] 延伸的频率向量,其中 f_max=1/(2*dt)。如果我们现在设计一个一维高斯滤波器(在傅立叶空间中)如下:
f_vec0=0;
sigma=1;
filter=exp( -(f_vec-f_vec0).^2./(2*sigma^2));
然后在傅立叶域中进行滤波:
f_signal=fftshift(fft(signal));
filt_signal=fftshift(ifft(f_signal.*filter));
因此,从我们应用的滤波器来看,sigma=1 意味着截止频率(我确定为滤波器最大值(即 1)的 1%)大约为 3 Hz:
cutoff_freq=f_vec(find(filter>=0.01,1,'last'))
将其转化为 2D 很简单,只需小心单位即可。对于图像来说,有像素作为位置单位,1/像素作为空间频率。 fspecial 函数生成预定义滤波器的二维矩阵。 fspecial的用法通常是这样的:
PSF = fspecial('gaussian',hsize,sigma);
Blurred = imfilter(Image,PSF,'symmetric','conv');
使用卷积就像在傅立叶域中进行乘法一样。傅里叶域中的 sigma 与位置域的 1/sigma 成正比,等等......