频域信号处理
快速傅里叶变换的应用非常广。用FFT(快速傅立叶变换)能将时域的数字信号转换为频域信号。转换为频域信号之后我们可以很方便地分析出信号的频率成分，在频域上进行处理，最终还可以将处理完毕的频域信号通过IFFT(逆变换)转换为时域信号，实现许多在时域无法完成的信号处理算法。
在python的scipy类库中提供了快速傅里叶变换包fftpack。下面的例子中低频的原始信号加入了高频噪声，我们可以通过快速傅里叶变换来分析其频谱，去掉噪声再重构信号。
import numpy as np from scipy import fftpack # The scipy.fftpack module allows to compute fast Fourier transforms import pylab as pl time_step = 0.02 # 信号采样间隔Δt period = 5. # 正弦信号周期T time_vec = np.arange(0, 20, time_step) # 生成采样点,采样次数N=20 # 模拟带高频噪声的信号 sig = np.sin(2 * np.pi / period * time_vec) + 0.5 * np.random.randn(time_vec.size) ''' fftfreq函数返回离散频率序列f, 频域间隔:Δf=1/(N*Δt) scipy.fftpack.fftfreq(n, d=1.0) f = [0, 1, ..., n/2-1, -n/2, ..., -1] / (d*n) if n is even f = [0, 1, ..., (n-1)/2, -(n-1)/2, ..., -1] / (d*n) if n is odd ''' sample_freq = fftpack.fftfreq(sig.size, d=time_step) sig_fft = fftpack.fft(sig) # compute the fast Fourier transform pidxs = np.where(sample_freq > 0) # only the positive part of the spectrum needs to be used freqs, power = sample_freq[pidxs], np.abs(sig_fft)[pidxs] freq = freqs[power.argmax()] # 去除高频噪声 sig_fft[np.abs(sample_freq) > freq] = 0 # 用傅里叶反变换，重构去除噪声的信号 main_sig = fftpack.ifft(sig_fft) pl.figure() pl.plot(freqs, power) pl.xlabel('Frequency [Hz]') pl.ylabel('plower') axes = pl.axes([0.3, 0.3, 0.5, 0.5]) pl.title('Peak frequency') pl.plot(freqs[:8], power[:8]) pl.setp(axes, yticks=[]) pl.figure() pl.plot(time_vec, sig) pl.plot(time_vec, main_sig, linewidth=3) pl.xlabel('Time [s]') pl.ylabel('Amplitude') pl.show()
下图可以看出在原始信号频率0.2Hz处，幅值最大，其它高于0.2Hz的都为噪声信号。用FFT方法消除噪声就是对含噪声信号的频谱进行处理，将噪声所在频段的$X(k)$值全部置零后，再对处理后的$X(k)$进行离散傅里叶逆变换(IFFT)可得原始信号的近似结果。这种方法要求噪声与信号的频谱不在同一频段，否则将很难处理。
下图中的粗线代表通过fft及ifft从带噪声的信号中提取的有用信号
- 二维离散傅里叶变换
一个图像为M×N的函数f(x,y)的离散傅里叶变换F(u,v)为:
$$F(u,v)=\frac{1}{MN}\sum_{x=0}^{M-1}\sum_{y=0}^{N-1}f(x,y)e^{-j2\pi (ux/M+vy/N)}$$
(u,v)=(0,0)位置的傅里叶变换为$F(0,0)=\frac{1}{MN}\sum_{x=0}^{M-1}\sum_{y=0}^{N-1}f(x,y)=\bar{f}(x,y)$，即f(x,y)的均值，原点(0,0)的傅里叶变换是图像的平均灰度。F(0,0)称为频率谱的直流分量，其它F(u,v)称为交流分量。傅里叶变换中出现的变量u和v通常称为频率变量，空间频率可以理解为等相位线在x,y坐标轴投影的截距的倒数。
相应的空间频率分别为：$u=\frac{1}{X},\quad v=\frac{1}{Y}$，对图像信号而言，空间频率是指单位长度内亮度作周期性变化的次数。
为了更好理解图像的傅里叶变换用下面几幅图来进行说明：
下图右边是原始图像，左边是通过傅里叶变换得到的图像（没有移频，即原点在图像的左上角）。由于原始图像的灰度只在Y轴上发生周期性变化，因此对应的傅里叶图像中只在左侧边缘出现一个亮点，代表原始图像Y轴上灰度变化的频率。如果右侧图像水平黑白条纹越密集（变化越快），则左侧图像上的亮点越远离原点向下；若右侧水平黑白条纹越稀疏（变化越慢），则左侧图像上的亮点越接近原点。
下面这幅图右侧的原始图像在X轴和Y轴上都有灰度的变化，因此对应的傅里叶图像上的亮点不再处于坐标轴上（边缘上）：
如果原始图像中存在两种频率，则显然其傅里叶变换的图像中会出现两个对应的离散点：
- 图像傅里叶变换的理解
把图像想象成沿着两个方向采集的信号。所以对图像同时进行X方向和Y方向的傅里叶变换，我们就会得到这幅图像的频域表示（频谱图），也就是图像梯度的分布图，当然频谱图上的各点与图像上各点并不存在一一对应的关系。二维傅里叶变换的坐标轴意义是频率，越靠近原点，频率越低，对应于图像中像素值变化速度比较慢的部分，一般是物体的主体、背景等。越远离原点，频率越高，对应于图像中像素值变化速度快的那部分，例如物体的边界。对一般图像来说靠近原点周围比较亮，远离原点的地方比较暗，也就是图像里低频部分分量多，高频分量少。一张图片中显然是边缘部分少，而大部分都是颜色相近/灰度相近的区域。
不同频率信息在图像结构中有不同的作用。图像的主要成分是低频信息，它形成了图像的基本灰度等级，对图像结构的决定作用较小；中频信息决定了图像的基本结构，形成了图像的主要边缘结构；高频信息形成了图像的边缘和细节，是在中频信息上对图像内容的进一步强化。傅立叶变换提供另外一个角度来观察图像，可以将图像从灰度分布转化到频率分布上来观察图像的特征。
- 利用傅里叶变换进行滤波
Numpy的FFT包可以帮助我们实现快速傅里叶变换。函数np.fft.fft2()可以对信号进行频率转换，输出结果是一个复杂的数组。函数的第一个参数是输入图像，要求是灰度格式。第二个参数是可选的，决定输出数组的大小。默认输出数组的大小和输入图像大小一样；如果输出结果比输入图像大，输入图像就需要在进行FFT前补0；如果输出结果比输入图像小的话，输入图像就会被切割。
import cv2 import numpy as np from matplotlib import pyplot as plt img = cv2.imread('input.png',0) # fft2() computes the n-dimensional discrete Fourier Transform f = np.fft.fft2(img) # Shift the zero-frequency component to the center of the spectrum fshift = np.fft.fftshift(f) magnitude_spectrum = 20*np.log(np.abs(fshift)) rows, cols = img.shape crow, ccol = rows/2 , cols/2 fshift[crow-40:crow+40, ccol-40:ccol+40] = 0 # remove the low frequencies by masking with a rectangular window of size 80x80 f_ishift = np.fft.ifftshift(fshift) img_back = np.fft.ifft2(f_ishift) img_back = np.abs(img_back) plt.subplot(131),plt.imshow(img, cmap = 'gray') plt.title('Input Image'), plt.xticks([]), plt.yticks([]) plt.subplot(132),plt.imshow(magnitude_spectrum, cmap = 'gray') plt.title('Magnitude Spectrum'), plt.xticks([]), plt.yticks([]) plt.subplot(133),plt.imshow(img_back, cmap = 'gray') plt.title('Image after HPF'), plt.xticks([]), plt.yticks([]) plt.show()
上图的结果显示高通滤波其实是一种边界检测操作。
OpenCV中相应的函数是cv2.dft()和cv2.idft()，和前面输出的结果一样，但是是双通道的。第一个通道是结果的实数部分，第二个通道是结果的虚数部分，输入图像要首先转换成np.float32格式。OpenCV中的函数cv2.dft() 和cv2.idft()要比Numpy快，但是Numpy函数更加用户友好。
当数组的大小为某些值时DFT的性能会更好。当数组的大小是2的幂时DFT效率最高；当数组的大小是 2，3，5的倍数时效率也会很高。所以如果你想提高代码的运行效率时，你可以修改输入图像的大小（补0）。对于OpenCV必须自己手动补0，但是Numpy只需要指定FFT运算的大小，它会自动补 0。那我们怎样确定最佳大小呢？ OpenCV提供了一个函数:cv2.getOptimalDFTSize()。它可以同时被cv2.dft()和np.fft.fft2()使用。
让我们使用 IPython中的魔法命令%timeit 来测试一下：
数组的大小从(186,295)变成了(192,300)，现在我们为它补 0，然后看看性能有没有提升。
现在我们看看Numpy的表现：
速度提高了3倍。我们再看看OpenCV的表现：
我们也会发现OpenCV的速度比Numpy还要快2倍。
在前面的部分我们实现了一个HPF（高通滤波），现在我们来实现LPF（低通滤波）将高频部分去除。其实就是对图像进行模糊操作。首先我们需要构建一个掩模，与低频区域对应的地方设置为 1, 与高频区域对应的地方设置为 0。
import cv2 import numpy as np from matplotlib import pyplot as plt img = cv2.imread('input.png',0) rows, cols = img.shape nrows = cv2.getOptimalDFTSize(rows) # Returns the optimal DFT size for a given vector size ncols = cv2.getOptimalDFTSize(cols) nimg = np.zeros((nrows,ncols)) nimg[:rows,:cols] = img dft = cv2.dft(np.float32(nimg), flags = cv2.DFT_COMPLEX_OUTPUT) dft_shift = np.fft.fftshift(dft) crow, ccol = nrows/2 , ncols/2 # create a mask first, center square is 1, remaining all zeros mask = np.zeros((nrows, ncols, 2), np.uint8) mask[crow-30:crow+30, ccol-30:ccol+30] = 1 # apply mask and inverse DFT fshift = dft_shift * mask f_ishift = np.fft.ifftshift(fshift) # Computes the inverse Discrete Fourier Transform of a 1D or 2D array img_back = cv2.idft(f_ishift) # You can also use cv2.cartToPolar() which returns both magnitude and phase img_back = cv2.magnitude(img_back[:,:,0],img_back[:,:,1]).reshape((nrows,ncols)) img_back = img_back[:rows,:cols] plt.subplot(121),plt.imshow(img, cmap = 'gray') plt.title('Input Image'), plt.xticks([]), plt.yticks([]) plt.subplot(122),plt.imshow(img_back, cmap = 'gray') plt.title('Image after LPF'), plt.xticks([]), plt.yticks([]) plt.show()
下面的一个例子利用二维傅里叶变换来消除图像上存在的周期性噪声。周期噪声一般产生于图像采集过程中的电气或电机的干扰，表现为图像中周期性的冲击，如下图：
通过观察傅立叶变换后的频谱图，我们就可以看出图像的能量分布，如果频谱图中暗的点数更多，那么实际图像是比较柔和的（因为各点与邻域差异都不大，梯度相对较小）；反之，如果频谱图中亮的点数多，那么实际图像一定是尖锐的，边界分明且边界两边像素差异较大的。对频谱移频到图像中心以后，可以看出图像的频率分布是以中心为圆心，对称分布的。将频谱移频到中心除了可以清晰地看出图像频率分布以外，还有一个好处，它可以分离出有周期性规律的干扰信号。比如正弦干扰，一幅带有正弦干扰，移频到中心的频谱图上可以看出除了中心以外还存在以某一点为中心，对称分布的亮点集合，这个集合就是干扰噪音产生的，这时可以很直观的通过在该位置放置带阻滤波器消除干扰。
若原点设在中心，其频谱能量集中分布在中心附近；若原点设在左上角，那么图像信号能量将集中在四个角上。这是由二维傅立叶变换本身性质决定的。变换之后的图像在原点平移之前四角是低频，最亮，平移之后中间部分是低频，最亮，亮度大说明低频的能量大。
这张登月照片上面有着很有规律的噪点，那么其FFT频谱图上面就会有非常规则的点。这些点就是噪声在频域空间的对应。
我们需要保留频谱图像中有用的部分，去掉高频噪声。下面的Python代码中设定了一个比例，保留这些有用信息，将高频区域的频率设为0，消除噪声影响。注意程序中没有进行移频操作。
import cv2 import numpy as np from matplotlib import pyplot as plt img = cv2.imread('moonlanding.png',0) rows, cols = img.shape nrows = cv2.getOptimalDFTSize(rows) # Returns the optimal DFT size for a given vector size ncols = cv2.getOptimalDFTSize(cols) f = np.fft.fft2(img, [nrows,ncols]) magnitude_spectrum_1 = 20 * np.log(1 + np.abs(f)) # Define the fraction of coefficients (in each direction) we keep keep_fraction = 0.1 # Set to zero all rows with indices between row*keep_fraction and # row*(1-keep_fraction): f[nrows*keep_fraction:nrows*(1-keep_fraction)] = 0 # Similarly with the columns: f[:, ncols*keep_fraction:ncols*(1-keep_fraction)] = 0 magnitude_spectrum_2 = 20 * np.log(1 + np.abs(f)) img_back = np.abs(np.fft.ifft2(f)) img_back = img_back[:rows,:cols] plt.subplot(221),plt.imshow(img, cmap = 'gray') plt.title('Input Image'), plt.xticks([]), plt.yticks([]) plt.subplot(222),plt.imshow(magnitude_spectrum_1, cmap = 'gray') plt.title('Magnitude Spectrum'), plt.xticks([]), plt.yticks([]) plt.subplot(223),plt.imshow(img_back, cmap = 'gray') plt.title('Reconstructed Image'), plt.xticks([]), plt.yticks([]) plt.subplot(224),plt.imshow(magnitude_spectrum_2, cmap = 'gray') plt.title('Filtered Spectrum'), plt.xticks([]), plt.yticks([]) plt.show()
在Mathematica软件中可以很方便的用函数ImagePeriodogram[image]来显示图像的频谱(shows the squared magnitude of the discrete Fourier transform (power spectrum) of image)。在该函数中使用可选参数Alignment可以设置零频率的位置，默认值是Alignment→Center，即图像中心频率最低。如果设置为Alignment→{Left,Top}，那么零频分量将被放置在图像左上角。
参考：
http://old.sebug.net/paper/books/scipydoc/frequency_process.html
Scipy : high-level scientific computing
How to use 2D Fourier analysis to clean the noise in an image
