03.fourier transform（傅立叶变换）

一：演示

原图

傅里叶变换

原本四个角比较亮(低频)，经过四个角平移到中心，中心会比较亮

中心亮的是低频；
四周暗的是高频；

将图像四周置0（暗的区域）后，傅里叶反变换图像：

将图像中心置0（亮的区域）后，傅里叶反变换图像：

二：什么是频域？

红色图像：时域图像；
蓝色图像：频域图像；

三：傅里叶变换的由来

三：傅里叶变换算法（DFT / FFT）

四：源码

fft2d.py

# -*- coding: utf-8 -*-
import numpy as np
import matplotlib.pyplot as pltdef rawFFT(Img, Wr, axis):if Img.shape[axis] == 2:pic = np.zeros(Img.shape, dtype=complex)if axis == 1:pic[:, 0] = Img[:, 0] + Img[:, 1] * Wr[:, 0]pic[:, 1] = Img[:, 0] - Img[:, 1] * Wr[:, 0]elif axis == 0:pic[0, :] = Img[0, :] + Img[1, :] * Wr[0, :]pic[1, :] = Img[0, :] - Img[1, :] * Wr[0, :]return picelse:pic = np.empty(Img.shape, dtype=complex)if axis == 1:A = rawFFT(Img[:, ::2], Wr[:, ::2], axis)B = rawFFT(Img[:, 1::2], Wr[:, ::2], axis)pic[:, 0:Img.shape[1] // 2] = A + Wr * Bpic[:, Img.shape[1] // 2:Img.shape[1]] = A - Wr * Belif axis == 0:A = rawFFT(Img[::2, :], Wr[::2, :], axis)B = rawFFT(Img[1::2, :], Wr[::2, :], axis)pic[0:Img.shape[0] // 2, :] = A + Wr * Bpic[Img.shape[0] // 2:Img.shape[0], :] = A - Wr * Breturn picdef FFT_1d(Img, axis):Wr = np.zeros(Img.shape, dtype=complex)if axis == 0:Wr = np.zeros((Img.shape[0] // 2, Img.shape[1]), dtype=complex)temp = np.array([np.cos(2 * np.pi * i / Img.shape[0]) - 1j * np.sin(2 * np.pi * i / Img.shape[0]) for i inrange(Img.shape[0] // 2)])for i in range(Wr.shape[1]):Wr[:, i] = tempelif axis == 1:Wr = np.zeros((Img.shape[0], Img.shape[1] // 2), dtype=complex)temp = np.array([np.cos(2 * np.pi * i / Img.shape[1]) - 1j * np.sin(2 * np.pi * i / Img.shape[1]) for i inrange(Img.shape[1] // 2)])for i in range(Wr.shape[0]):Wr[i, :] = tempreturn rawFFT(Img, Wr, axis)def iFFT_1d(Img, axis):Wr = np.zeros(Img.shape, dtype=complex)if axis == 0:Wr = np.zeros((Img.shape[0] // 2, Img.shape[1]), dtype=complex)temp = np.array([np.cos(2 * np.pi * i / Img.shape[0]) + 1j * np.sin(2 * np.pi * i / Img.shape[0]) for i inrange(Img.shape[0] // 2)])for i in range(Wr.shape[1]):Wr[:, i] = tempelif axis == 1:Wr = np.zeros((Img.shape[0], Img.shape[1] // 2), dtype=complex)temp = np.array([np.cos(2 * np.pi * i / Img.shape[1]) + 1j * np.sin(2 * np.pi * i / Img.shape[1]) for i inrange(Img.shape[1] // 2)])for i in range(Wr.shape[0]):Wr[i, :] = tempreturn rawFFT(Img, Wr, axis) * (1.0 / Img.shape[axis])def FFT_2d(Img):'''only for gray scale 2d-img. otherwise return 0 img with the same size of Img:param Img: img to be fourior transform:return: img been transformed'''imgsize = Img.shapepic = np.zeros(imgsize, dtype=complex)if len(imgsize) == 2:N = 2while N < imgsize[0]:N = N << 1num1 = N - imgsize[0]N = 2while N < imgsize[1]:N = N << 1num2 = N - imgsize[1]pic = FFT_1d(np.pad(Img, ((num1 // 2, num1 - num1 // 2), (0, 0)), 'edge'), 0)[num1 // 2:num1 // 2 + imgsize[0], :]pic = FFT_1d(np.pad(pic, ((0, 0), (num2 // 2, num2 - num2 // 2)), 'edge'), 1)[:,num2 // 2:num2 // 2 + imgsize[1]]return picdef iFFT_2d(Img):'''only for gray scale 2d-img. otherwise return 0 img with the same size of Img:param Img: img to be fourior transform:return: img been transformed'''imgsize = Img.shapepic = np.zeros(imgsize, dtype=complex)if len(imgsize) == 2:N = 2while N < imgsize[0]:N = N << 1num1 = N - imgsize[0]N = 2while N < imgsize[1]:N = N << 1num2 = N - imgsize[1]pic = iFFT_1d(np.pad(Img, ((num1 // 2, num1 - num1 // 2), (0, 0)), 'edge'), 0)[num1 // 2:num1 // 2 + imgsize[0],:]  # ,constant_values=(255,255)pic = iFFT_1d(np.pad(pic, ((0, 0), (num2 // 2, num2 - num2 // 2)), 'edge'), 1)[:,num2 // 2:num2 // 2 + imgsize[1]]  # ,constant_values=(255,255)return picif __name__ == "__main__":img = plt.imread('./img/1.jpg')img = img.mean(2)  # grayplt.imshow(img.astype(np.uint8), cmap='gray')plt.axis('off')plt.show()F1 = np.fft.fft2(img[:256, :256])F2 = FFT_2d(img[:256, :256])print((abs(F1 - F2) < 0.0000001).all())F1 = np.fft.ifft2(F1)F2 = iFFT_2d(F2)print((abs(F1 - F2) < 0.0000001).all())F2 = np.abs(F2)F2[F2 > 255] = 255plt.imshow(F2.astype(np.uint8), cmap='gray')plt.axis('off')plt.show()F1 = np.abs(F1)F1[F1 > 255] = 255plt.imshow(F1.astype(np.uint8), cmap='gray')plt.axis('off')plt.show()

Presentation.py

import matplotlib.pyplot as plt
import numpy as npif __name__ == '__main__':img = plt.imread('./img/raccoon.jpg')plt.imshow(img, cmap=plt.cm.gray)plt.axis('off')plt.show()img = img.mean(axis=-1)# plt.imsave('gray_raccoon.jpg', np.dstack((img.astype(np.uint8), img.astype(np.uint8), img.astype(np.uint8))))img = np.fft.fft2(img)img = np.fft.fftshift(img)fourier = np.abs(img)magnitude_spectrum = np.log(fourier)plt.imshow(magnitude_spectrum.astype(np.uint8), cmap=plt.cm.gray)plt.axis('off')plt.show()  # image after fourier transform# plt.imsave('fourier_raccoon.jpg', 14*np.dstack((magnitude_spectrum.astype(np.uint8),magnitude_spectrum.astype(np.uint8),magnitude_spectrum.astype(np.uint8))))x, y = img.shapelowF = np.zeros((x, y))lowF = lowF.astype(np.complex128)window_shape = (20, 20)lowF[int(x / 2) - window_shape[0]:int(x / 2) + window_shape[0],int(y / 2) - window_shape[1]:int(y / 2) + window_shape[1]] = \img[int(x / 2) - window_shape[0]:int(x / 2) + window_shape[0],int(y / 2) - window_shape[1]:int(y / 2) + window_shape[1]]lowF_im = np.fft.ifft2(lowF)lowF_im = np.abs(lowF_im)lowF_im[lowF_im > 255] = 255plt.imshow(lowF_im.astype(np.uint8), cmap='gray')plt.axis('off')plt.show()# plt.imsave('LowF_raccoon.jpg', np.dstack((lowF_im.astype(np.uint8), lowF_im.astype(np.uint8), lowF_im.astype(np.uint8))))highF = np.zeros((x, y))highF = highF.astype(np.complex128)window_shape = (370, 370)highF[0:window_shape[0], :] = img[0:window_shape[0], :]highF[x - window_shape[0]:x, :] = img[x - window_shape[0]:x, :]highF[:, 0:window_shape[1]] = img[:, 0:window_shape[1]]highF[:, y - window_shape[1]:y] = img[:, y - window_shape[1]:y]highF_im = np.fft.ifft2(highF)highF_im = np.abs(highF_im)highF_im[highF_im > 255] = 255plt.imshow(highF_im.astype(np.uint8), cmap='gray')plt.axis('off')plt.show()# plt.imsave('HighF_raccoon.jpg', np.dstack((highF_im.astype(np.uint8), highF_im.astype(np.uint8), highF_im.astype(np.uint8))))

LiveFourierTransform.py

import numpy as np
import time
import cv2def LiveFFT(rows, columns, numShots=1e5, nContrIms=30):imMin = .004  # minimum allowed value of any pixel of the captured imagecontrast = np.concatenate((np.zeros((10, 1)), np.ones((10, 1))), axis=1)  # internal use.window_name = 'Livefft by lixin'vc = cv2.VideoCapture(0)  # camera devicecv2.namedWindow(window_name, 0)  # 0 makes it work a bit bettercv2.resizeWindow(window_name, 1024, 768)  # this doesn't keeprval, frame = vc.read()# we need to wait a bit before we get decent imagesprint("warming up camera... (.1s)")time.sleep(.1)rval = vc.grab()rval, frame = vc.retrieve()# determine if we are not asking too muchframeShape = np.shape(frame)if rows > frameShape[1]:rows = frameShape[1]if columns > frameShape[0]:columns = frameShape[0]# calculate cropvCrop = np.array([np.ceil(frameShape[0] / 2. - columns / 2.),np.floor(frameShape[0] / 2. + columns / 2.)], dtype=int)hCrop = np.array([np.ceil(frameShape[1] / 2. - rows / 2.),np.floor(frameShape[1] / 2. + rows / 2.)], dtype=int)# start image cleanup with something like this:# for a running contrast of nContrIms framescontrast = np.concatenate((np.zeros((nContrIms, 1)),np.ones((nContrIms, 1))),axis=1)Nr = 0# main loopwhile Nr <= numShots:a = time.time()Nr += 1contrast = work_func(vCrop, hCrop, vc, imMin, window_name, contrast)print('framerate = {} fps \r'.format(1. / (time.time() - a)))# stop cameravc.release()def work_func(vCrop, hCrop, vc, imMin, figid, contrast):# read imagerval = vc.grab()rval, im = vc.retrieve()im = np.array(im, dtype=float)# crop imageim = im[vCrop[0]: vCrop[1], hCrop[0]: hCrop[1], :]# pyramid downscaling# im = cv2.pyrDown(im)# reduce dimensionalityim = np.mean(im, axis=2, dtype=float)# make sure we have no zerosim = (im - im.min()) / (im.max() - im.min())im = np.maximum(im, imMin)Intensity = np.abs(np.fft.fftshift(np.fft.fft2(im))) ** 2Intensity += imMin# kill the center lines for higher dynamic range# by copying the next row/column# h, w = np.shape(Intensity)# Intensity[(h / 2 - 1):(h / 2 + 1), :] = Intensity[(h / 2 + 1):(h / 2 + 3), :]# Intensity[:, (w / 2 - 1):(w / 2 + 1)] = Intensity[:, (w / 2 + 1):(w / 2 + 3)]# running average of contrast##circshift contrast matrix upcontrast = contrast[np.arange(1, np.size(contrast, 0) + 1) % np.size(contrast, 0), :]##replace bottom values with new values for minimum and maximumcontrast[-1, :] = [np.min(Intensity), np.max(Intensity)]maxContrast = 1minContrast = 7  # to be modify# openCV drawvmin = np.log(contrast[:, 0].mean()) + minContrastvmax = np.log(contrast[:, 1].mean()) - maxContrastIntensity = (np.log(Intensity + imMin) - vmin) / (vmax - vmin)Intensity = Intensity.clip(0., 1.)# Intensity = (Intensity - Intensity.min()) / (Intensity.max() - Intensity.min())time.sleep(.01)cv2.imshow(figid, np.concatenate((im, Intensity), axis=1))cv2.waitKey(1)return contrastif __name__ == '__main__':LiveFFT(400, 400)

摘录：会飞的吴克笔记
b站有视频！

03.fourier transform（傅立叶变换）相关推荐

使用 scipy.fft 进行Fourier Transform：Python 信号处理
摘要:Fourier transform 是一个强大的概念,用于各种领域,从纯数学到音频工程甚至金融. 本文分享自华为云社区<使用 scipy.fft 进行Fourier Transform:P ...
深度学习开源框架系列：基础算法之傅立叶变换：1：概要介绍
傅立叶变换时数字信号处理的重要方法之一,是法国数学家傅立叶在1807年在法国科学学会上发表的一篇文章中所提出的,在文章中使用了正弦函数描述温度分布,而且提出了一个著名的论断:任何连续性的周期信号都可以 ...
傅立叶变换(Fourier Transform)分析理解
引言关于傅立叶变换,无论是书本还是在网上可以很容易找到关于傅立叶变换的描述,但是大都是些故弄玄虚的文章,太过抽象,尽是一些让人看了就望而生畏的公式的罗列,让人很难能够从感性上得到理解,最近,我偶尔从 ...
matlab prefourier,用matlab求单位阶跃函数的傅立叶变换fourier变换
傅立叶变换常常应用于在信号处理,其在信号处理中有着很重要的位置.下面我们就介绍一种最简单的傅立叶变换,即单位阶跃函数的傅立叶变换.然后用matlab实现单位阶跃函数的傅立叶变换. 1.傅立叶变换简介 ...
3D数学之快速傅立叶变换(Fast Fourier Transform-FFT)
最近写毕业论文,想做一个基于FFT的海洋,但是时间比较紧,不知道来不来得及去实现它,但是FFT从我研一开始就听说过了,一直没去特意去研究过它,趁此机会,我想在此把这个知识点搞定! 这篇文章可以耗时比较 ...
OpenCV 【十七】离散傅立叶变换
目录 1 key 2 原理 3 实例 3代码 4运行结果 5应用举例 1 key 什么是傅立叶变换及其应用? 如何使用OpenCV提供的傅立叶变换? 相关函数的使用,如: copyMakeBorder ...
理解离散傅立叶变换（一）
理解离散傅立叶变换(一) ------傅立叶变换的由来关于傅立叶变换,无论是书本还是在网上可以很容易找到关于傅立叶变换的描述,但是大都是些故弄玄虚的文章,太过抽象,尽是一些让人看了就望而生畏的公式 ...
UA OPTI512R 傅立叶光学导论17 离散傅立叶变换简介
UA OPTI512R 傅立叶光学导论17 离散傅立叶变换简介 DFT及其矩阵形式 DFT的性质上一讲提到对连续波形f(x)f(x)f(x)做周期性采样时可以用采样函数来表示采样结果: fS(x)= ...
为什么要进行傅立叶变换？傅立叶变换究竟有何意义？如何用Matlab实现快速傅立叶变换？
https://www.douban.com/note/164400821/ 写在最前面:本文是我阅读了多篇相关文章后对它们进行分析重组整合而得,绝大部分内容非我所原创.在此向多位原创作者致敬!!! ...