信号处理在频域中的应用
快速傅里叶变换的应用非常广。用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
SciPy Lecture Notes 中文版
图像傅里叶变换
傅里叶变换在图像处理中的应用
知乎-傅里叶变换有哪些具体的应用?
How to use 2D Fourier analysis to clean the noise in an image
An Intuitive Explanation of Fourier Theory
What does frequency domain denote in case of images?
上一篇: [Python图像处理] 32. 详细总结傅里叶变换去噪与霍夫变换特征识别(附带万字解析)
下一篇: 数字图像处理中的图像变换:二维离散傅里叶正逆变换(FFT2/IFFT2)、离散余弦正逆变换(DCT2/IDCT2)、频谱正逆平移(FFTSHIFT/IFFTSHIFT)、幅度谱与相位谱的例题与分析
推荐阅读
-
Java 8新特性探究(十三)JavaFX 8新特性以及开发2048游戏-JavaFX历史## 跟java在服务器端和web端成绩相比,桌面一直是java的软肋,于是Sun公司在2008年推出JavaFX,弥补桌面软件的缺陷,请看下图JavaFX一路走过来的改进 从上图看出,一开始推出时候,开发者需使用一种名为JavaFX Script的静态的、声明式的编程语言来开发JavaFX应用程序。因为JavaFX Script将会被编译为Java bytecode,程序员可以使用Java代码代替。 JavaFX 2.0之后的版本摒弃了JavaFX Script语言,而作为一个Java API来使用。因此使用JavaFX平台实现的应用程序将直接通过标准Java代码来实现。 JavaFX 2.0 包含非常丰富的 UI 控件、图形和多媒体特性用于简化可视化应用的开发,WebView可直接在应用中嵌入网页;另外 2.0 版本允许使用 FXML 进行 UI 定义,这是一个脚本化基于 XML 的标识语言。 从JDK 7u6开始,JavaFx就与JDK捆绑在一起了,JavaFX团队称,下一个版本将是8.0,目前所有的工作都已经围绕8.0库进行。这是因为JavaFX将捆绑在Java 8中,因此该团队决定跳过几个版本号,迎头赶上Java 8。 ##JavaFx8的新特性 ## ###全新现代主题:Modena 新的Modena主题来替换原来的Caspian主题。不过在Application的start方法中,可以通过setUserAgentStylesheet(STYLESHEET_CASPIAN)来继续使用Caspian主题。 参考http://fxexperience.com/2013/03/modena-theme-update/ ###JavaFX 3D 在JavaFX8中提供了3D图像处理API,包括Shape3D (Box, Cylinder, MeshView, Sphere子类),SubScene, Material, PickResult, LightBase (AmbientLight 和PointLight子类),SceneAntialiasing等。Camera类也得到了更新。从JavaDoc中可以找到更多信息。 ###富文本 强化了富文本的支持 ###TreeTableView ###日期控件DatePicker 增加日期控件 ###用于 CSS 结构的公共 API
-
信号处理在频域中的应用
-
基于FNN网络的电影评论情感分析在自然语言处理(NLP)中的应用
-
基于预训练模型的机器阅读理解在自然语言处理(NLP)中的应用
-
使用PaddleHub的文本审核服务在自然语言处理中的应用
-
深入理解JSP的九大内置对象及其在四个作用域中的应用
-
小白也能懂!广义特征值分解在图像处理中的应用教程
-
神经网络在多元回归数据分析中的应用与数据处理技巧
-
实战探索:Apache Hudi在流处理与批处理场景中的应用实例
-
【摩尔线程+Colossal-AI强强联手】MusaBert登上CLUE榜单TOP10:技术细节揭秘 - 技术实力:摩尔线程凭借"软硬兼备"的技术底蕴,让MusaBert得以从底层优化到顶层。其内置多功能GPU配备AI加速和并行计算模块,提供了全面的AI与科学计算支持,为AI推理和低资源条件下的大模型训练等场景带来了高效、经济且环保的算力。 - 算法层面亮点:依托Colossal-AI AI大模型开发系统,MusaBert在训练过程中展现出了卓越的并行性能与易用性,特别在预处理阶段对DataLoader进行了优化,适应低资源环境高效处理海量数据。同时,通过精细的建模优化、领域内数据增强以及Adan优化器等手段,挖掘和展示了预训练语言模型出色的语义理解潜力。基于MusaBert,摩尔线程自主研发的MusaSim通过对比学习方法微调,结合百万对标注数据,MusaSim在多个任务如语义相似度、意图识别和情绪分析中均表现出色。 - 数据资源丰富:MusaBert除了自家高质量语义相似数据外,还融合了悟道开源200GB数据、CLUE社区80GB数据,以及浪潮公司提供的1TB高质量数据,保证模型即便在较小规模下仍具备良好性能。 当前,MusaBert已成功应用于摩尔线程的智能客服与数字人项目,并广泛服务于语义相似度、情绪识别、阅读理解与声韵识别等领域。为了降低大模型开发和应用难度,MusaBert及其相关高质量模型代码已在Colossal-AI仓库开源,可快速训练优质中文BERT模型。同时,通过摩尔线程与潞晨科技的深度合作,仅需一张多功能GPU单卡便能高效训练MusaBert或更大规模的GPT2模型,显著降低预训练成本,进一步推动双方在低资源大模型训练领域的共享目标。 MusaBert荣登CLUE榜单TOP10,象征着摩尔线程与潞晨科技联合研发团队在中文预训练研究领域的领先地位。展望未来,双方将携手探索更大规模的自然语言模型研究,充分运用上游数据资源,产出更为强大的模型并开源。持续强化在摩尔线程多功能GPU上的大模型训练能力,特别是在消费级显卡等低资源环境下,致力于降低使用大模型训练的门槛与成本,推动人工智能更加普惠。而潞晨科技作为重要合作伙伴,将继续发挥关键作用。