1. 傅里叶变换原理

在数学中进行傅里叶变换为连续模拟信号，通常来说：
二维连续函数f(x,y)的傅里叶正变换为：

相应的傅里叶逆变换公式为：

但是在计算机领域，计算机一般处理的是数字信号，只能进行有限次计算，因此将这种受限下的傅里叶变换成为离散傅里叶变换（Discrete Fourier Transform,DFT)。

二维离散函数f(x,y)的傅里叶正变换的公式如下：

相应的傅里叶逆变换的公式如下：

2. 傅里叶变换性质

这里列出来书本上介绍的关于傅里叶变换的一些性质，仅供参考，不做详述。

3. 傅里叶变换python实现

3.1 傅里叶正变换实现

参考上述傅里叶正变换给出的公式，这里给出其python实现如下：

def DFT2D(x, shift=True):'''Discrete space fourier transformx: Input matrix'''pi2 = 2*np.piN1, N2 = x.shapeX = np.zeros((N1, N2), dtype=np.complex64)n1, n2 = np.mgrid[0:N1, 0:N2]for w1 in range(N1):for w2 in range(N2):j2pi = np.zeros((N1, N2), dtype=np.complex64)j2pi.imag = pi2*(w1*n1/N1 + w2*n2/N2)X[w1, w2] = np.sum(x*np.exp(-j2pi))if shift:X = np.roll(X, N1//2, axis=0)X = np.roll(X, N2//2, axis=1)return X

3.2 傅里叶逆变换实现

参考上述傅里叶逆变换给出的公式，这里给出其python实现如下：

def iDFT2D(X, shift=True):'''Inverse discrete space fourier transformX: Complex matrix'''pi2 = 2*np.piN1, N2 = X.shapex = np.zeros((N1, N2))k1, k2 = np.mgrid[0:N1, 0:N2]if shift:X = np.roll(X, -N1//2, axis=0)X = np.roll(X, -N2//2, axis=1)for n1 in range(N1):for n2 in range(N2):j2pi = np.zeros((N1, N2), dtype=np.complex64)j2pi.imag = pi2*(n1*k1/N1 + n2*k2/N2)x[n1, n2] = abs(np.sum(X*np.exp(j2pi)))return 1/(N1*N2)*x

3.3 测试用例

这里给出测试用例，如下：

import matplotlib.pyplot as plt
import numpy as np
import imageio
from _utils import *
image = imageio.imread('./sample/cameraman.png')
s = 4
image = image[::s, ::s]/255
N1, N2 = image.shape
IMAGE = DFT2D(image)
xX = np.array([image, np.log10(1 + abs(IMAGE))])
panel(xX, [2, 1], text_color='green',texts=['Input image', 'Spectrum'])
image_ = iDFT2D(IMAGE)
Xx_ = np.array([np.log10(1 + abs(IMAGE)), image_])
panel(Xx_, [2, 1], text_color='green',texts=['Spectrum', 'Reconstructed image'])

傅里叶正变换运行结果如下：

傅里逆正变换运行结果如下：

4. 傅里叶变换cuda实现

上述Python代码的运行时间复杂度为 O(N^2)，图像较大时，运行贼慢。这里考虑用cuda对其加速，搜索了半天，发现有大佬曾经做过类似的实现，并给出了源码，这里直接参考其实现对其进行简单的封装，可以调用完成图像的傅里叶变换以及视频的傅里叶变换。

这里给出FFT.cu的核心代码片段：

#include "fft.h"
#include "cuda_runtime.h"__device__ unsigned char getr(float x) {return (tanh((x - 0.375f) * 6.0f) + 1.0f) * 127.0f;
}
__device__ unsigned char getg(float x) {return (tanh((x - 0.6250f) * 6.0f) + 1.0f) * 127.0f;
}
__device__ unsigned char getb(float x) {return (exp(-20.0f * (x - 0.25f) * (x - 0.25f) - 2.0f * exp(-(x + 0.05f) * (x + 0.05f) * 144.0f)) * 0.5f + 1.0f + tanh((x - 0.875f) * 6.0f)) * 127.0f;
}
__global__ void imgfill(float2* d_k, uchar3* d_img,int size)
{int x = threadIdx.x + blockIdx.x * blockDim.x;int y = threadIdx.y + blockIdx.y * blockDim.y;int imgx, imgy;imgx = (x >= size / 2) ? x - size / 2 : x + size / 2;imgy = (y >= size / 2) ? y - size / 2 : y + size / 2;float2 k = d_k[y * size + x];float in = k.x * k.x + k.y * k.y;in = log(in * (1.0f / 256.0f/size) + 0.8f) * 0.07f;uchar3 c;c.x = getb(in);c.y = getg(in);c.z = getr(in);d_img[imgy * size + imgx] = c;
}__global__ void fill(float2* d_x, uchar3* d_8uc3,int size,int w,int h) {int x = threadIdx.x + blockIdx.x * blockDim.x;int y = threadIdx.y + blockIdx.y * blockDim.y;int imgx, imgy;float cx, cy;unsigned char r;if (x >= size / 2 + w / 2) {imgx = 0;cx = size - x;cx = exp(-cx * cx * (1.0f / 1024.0f));}else if (x < size / 2 + w / 2 && x >= w) {imgx = w - 1;cx = x - w;cx = exp(-cx * cx * (1.0f / 1024.0f));}else {imgx = x;cx = 1.0f;}if (y >= size / 2 + h / 2) {imgy = 0;cy = size - y;cy = exp(-cy * cy * (1.0f / 1024.0f));}else if (y < size / 2 + h / 2 && y >= h) {imgy = h - 1;cy = y - h;cy = exp(-cy * cy * (1.0f / 1024.0f));}else {imgy = y;cy = 1.0f;}r = d_8uc3[imgy * w + imgx].x;d_x[y * size + x].x = r * cx * cy;d_x[y * size + x].y = 0;
}void fft_tranformer(uchar3 * d_8uc3,float2 * d_x,float2 * d_k,uchar3 *d_img,cufftHandle *fftPlan,unsigned char * pframe,unsigned  char * pDst,int width,int height,int size)
{cudaMemcpy(d_8uc3, pframe, width * height * 3, cudaMemcpyHostToDevice);fill << < dim3(size / 128, size, 1), dim3(128, 1, 1) >> > (d_x, d_8uc3,size,width,height);cufftExecC2C(*fftPlan, d_x, d_k, CUFFT_FORWARD);imgfill << < dim3(size / 128, size, 1), dim3(128, 1, 1) >> > (d_k, d_img,size);cudaMemcpy(pDst, d_img, size * size * 3, cudaMemcpyDeviceToHost);
}

5. cuda实现效果

5.1 静态效果

5.2 动态效果

5.3 完整带音乐视频效果

点我

6. 代码

完整代码可公众号内回复DFT 即可获取

7. 参考

链接一

链接二

关注公众号《AI算法之道》，获取更多AI算法资讯.

【深度好文】二维图像傅里叶变换 YYDS相关推荐

【深度好文】二维图像离散余弦变换
1.引言上节分享了二维图像离散傅里叶变换,本节来继续讲频域空间的另一种变换–二维离散余弦变换(Discrete Cosine Transform,DCT).从运算方式上来讲,离散傅里叶变换计算的对象 ...
第4章 Python 数字图像处理(DIP) - 频率域滤波5 - 二变量函数的傅里叶变换、图像中的混叠、二维离散傅里叶变换及其反变换
目录二变量函数的傅里叶变换二维冲激及其取样性质二维连续傅里叶变换对二维取样和二维取样定理图像中的混叠二维离散傅里叶变换及其反变换二变量函数的傅里叶变换二维冲激及其取样性质两个连续变量 ...
图像傅里叶变换（二维离散傅里叶变换）
图像傅里叶变换二维离散傅里叶变换是将图像从空间域转至频域,在图像增强.图像去噪.图像边缘检测.图像特征提取.图像压缩等等应用中都起着极其重要的作用.理论基础是任意函数都可以表示成正弦函数的线性组合的 ...
ICLR 2021|基于GAN的二维图像无监督三维形状重建
2D GAN知道3D形状吗?基于GAN的二维图像无监督三维形状重建论文.代码地址:在公众号「计算机视觉工坊」,后台回复「二维图像GAN」,即可直接下载. 摘要: 自然图像是三维物体在二维图像平面上的 ...
时间序列转二维图像方法及其应用研究综述
目录 1 前言 2 方法综述 2.1 时频分析法 2.1.1 短时傅里叶变换 2.1.2 小波变换 2.1.3 希尔伯特-黄变换 2.2 图像编码方法 2.2.1 格兰姆角场 2.2.2 马尔可夫转移 ...
【OpenCV 例程200篇】73. 二维连续傅里叶变换
[OpenCV 例程200篇]73. 二维连续傅里叶变换欢迎关注『OpenCV 例程200篇』系列,持续更新中欢迎关注『Python小白的OpenCV学习课』系列,持续更新中 2.1 二维 ...
程序员用「美貌」突破二维图像的人脸识别
GitChat 作者:于航原文: 如何利用"女装术"突破基于二维图像的人脸识别关注微信公众号:「GitChat 技术杂谈」一本正经的讲技术 [不要错过文末彩蛋] 首先声明,这 ...
C++实现二维离散傅里叶变换
在上一篇文章<C++实现一维离散傅里叶变换>中,我们介绍了一维信号傅立叶变换的公式和C++实现,并阐述了频域幅值的意义. 一维傅立叶变换只适用于一维信号,例如音频数据.心脑电图等. 在图像 ...
二维图像的傅立叶变换
摘要:二维图像的傅立叶变换,与一维傅立叶相比,在理解上要抽象很多.我在网上找了几篇相对较好的文章,并用matlab自己做了几个实验图像,希望能对大家理解二维图像的傅立叶变换有所帮助. 关键字:二维傅立 ...
AE插件-二维图像创建三维视觉错觉场景特效 Autostereogram
Autostereogram mac版带给大家,非常旨在从二维图像创建3D场景的视觉错觉,您只需要一张深度图和可能的纹理. AE中创建自动立体图-旨在从二维图像创建3D场景的视觉错觉,您只需要一张深度 ...

【深度好文】二维图像傅里叶变换 YYDS