c++ 提取傅里叶描述子_如何在Unity实现二维傅里叶变换与逆变换？

编者按 傅里叶变换是一种分析信号的方法，它可分析信号的成分，也可用这些成分合成信号。那么，你知道它也可以应用在游戏开发里吗？本文作者将分享一个具体工程上的完整案例，看看它是如何在Unity中实现的吧。

文 | lingzerg

(本文内容由知乎专栏“番茄的游戏开发世界”提供，转载请征得同意。文章仅为作者观点，不代表GWB立场)

在网上看到很多关于傅里叶变换的内容, 但是没找到具体工程上完整的一个例子。

例如把一个纹理转化为频谱图和相位，然后利用频谱和相位在转化回来。

于是就自己做一个好了。如果有不对之处请使劲喷。

然后如果你比较熟悉只想看工程部分的内容, 可以酌情跳过。

先看一眼结果:

放出Git 地址:

https://github.com/lingzerg/LingzergDemo

再看一眼流程:

然后我们慢慢解释。

什么是信号?

先把思路捋清楚。

讲傅里叶变换的定义有很多, 推荐大家各种搜一下。

比如3B1B的傅里叶变换:
https://www.bilibili.com/video/BV1pW411J7s8

比如这篇文章:
https://zhuanlan.zhihu.com/p/19763358

我简单描述下, 正常情况下，我们看到的信号是这样的:

这种信号你可以想象是声音，在时间轴上连续不断变换，自然界中的信号我们认为他是连续的，比如这样:

横轴t就是时间，所以这个叫做时域

这就是一个信号在时间轴上不断变化，所以这种信号我们叫一个信号的时域表示。

然后为了让计算机可以记录描述这种信号，我们会这样:

对信号进行离散采样

以一个固定间隔对信号进行采样，得到的数据包含这个信号在当前采样点的振幅和频率。

公式描述是:

f是频率(frequency), A是振幅(Amplitude), phase是相位,x是时间, 得到的是振幅。

所以我们只要是讨论计算机领域的信号, 就一定是离散的。

推荐大家去坐标系中手动拉一拉这些变量感受下这些变量的作用:

https://www.desmos.com/calculator/nyyk0mzlr6

什么是傅里叶变换?

我们拿到一个信号，有时候想过滤掉高频部分，例如通话降噪，还有一些不想要的低频部分。

这时我们可以用复数坐标系去描述一个信号, 让信号的另外一个特征暴露出来,以便于我们处理。

复数坐标系:

使用复数平面描述信号有什么好处呢?

在复数和实数构建的坐标系中，信号的时间周期就变成了一个环，信号采样就是不断围绕这个圈采样。

模长就是振幅，那问题就在于,我们如何把一个我们抽样得到的信号样本转换成复数形式?

先把欧拉公式扔这里: 

然后欧拉公式实际上是逆时针旋转的,而正变换是顺时针的。

所以我们需要改成

而逆变换就要换回来。

接着把频率  和时间  带入欧拉公式另一侧, 除以  是因为N次采样。

再乘以2π 转成弧度:

就如同公式里的描述:

而这个又可以展开写一篇了, 所以详细描述我推荐诸位看我上面发的那两个链接。

最后我们反正知道把一个信号的采样结果带入离散傅里叶变换的公式就可以得到一组复数的描述。

在这个公式里，实际上 x 就是时间点, u就是频率,小  就是每个单位时间的采样点采样得到的振幅,最后会得到一个复数, 而逆变换就是反过来:

大  就是傅里叶变换得到的复数。

二维傅里叶变换?

一维傅里叶变换和二维的差别，就在于二维还要在套一层，你需要遍历整个二维空间。

二维的离散逆傅里叶变换:

而在图像处理中, 二维傅里叶变换上任意一个像素可以认为是一整个平面波:

还可以把图片上的每个像素看成一个在x轴上方的信号，比如这样:

一个sin(x)+1的信号

接着将每个采样点带入公式就可以得到每个像素对应的复数。

我开始做傅里叶变换的时候非常奇怪，我感觉傅里叶变换和逆变换的公式之间并没有频域和相位啊，那这个流程到底是什么?

我先上一张图然后慢慢解释。

这个基本上就是傅里叶变换的完整流程了。

也就是说, 傅里叶变换原始公式里是只有复数的,但是通过复数我们可以得到每个像素点的膜长, 还有弧度:

你不做这一步，直接把复数的实部和虚部分别保存在纹理的r和g通道也是完全没问题的。

但是这就失去了对二维图片进行傅里叶变换的意义。

我们只有获取了频域和时域,才能根据频域和时域处理图片。

我看过很多demo, 都是直接保存复数,而没有对频域和相位进行独立的转换。

那我们怎么在Unity里实现呢?

首先我选择用compute shader来实现这个。

因为感觉CS比较帅吧。

在Cshader中我们创建一个复数的结构体,并写好常用方法:

// Complexstruct Complex{    float real;    float imaginary;};Complex getComplex(float num1, float num2){    Complex c;    c.real = num1;    c.imaginary = num2;    return c;}Complex getAdd(Complex c, Complex o){    c.real = c.real + o.real;    c.imaginary = c.imaginary + o.imaginary;    return c;}Complex getSubtract(Complex c,Complex o){    c.real = c.real - o.real;    c.imaginary = c.imaginary - o.imaginary;    return c;}Complex getMultiply(Complex c,Complex o){    float tmp = c.real * o.real - c.imaginary * o.imaginary;    c.imaginary = c.real * o.imaginary + c.imaginary * o.real;    c.real = tmp;    return c;}Complex getMultiplyFloat(Complex c,float x){    c.real = c.real *x;    c.imaginary = c.imaginary * x;    return c;}//欧拉公式 e^(ix)=cosx+isinxComplex getEuler(float x) {    Complex c;    c.real = cos(x);    c.imaginary = sin(x);    return c;}int2 getCoordinate(int x, int y) {     x = x - 31;     y = y - 31;     return int2(x,y);}

然后我们修改一下CS的核心, 按照公式做积分,详细描述在注释中:

[numthreads(16,16,1)]void Fourier (uint3 id : SV_DispatchThreadID){    float flag = -1 * 2 * 3.14 / TextureSize; // 把e ^ 2π /M 放到这里, 没必要每次循环都算这个    int2 coord = id.xy;//getCoordinate(id.x,id.y); //替换这里可以把复平面移动到中间    Complex sumRow;//定义一个复数,用与积分    sumRow.real = 0;    sumRow.imaginary = 0;    for(int i = 0; i< TextureSize; i++) { // 循环当前行        Complex complexU = getEuler(flag*i*coord.x); //先算当前行对应列的首位像素的欧拉数, 下面列积分后要乘这个数        Complex sumColumn;// 定义一个复数,用于列的积分        sumColumn.real = 0;        sumColumn.imaginary = 0;        for(int j = 0; j < TextureSize; j++) { // 循环每行对应的列            Complex complexV = getEuler(flag*j*coord.y); // 计算当前像素的欧拉数 - 2 π * v * y / Len            complexV = getMultiplyFloat(complexV, originalImg[int2(i,j)].r); // 用欧拉数乘以当前的振幅, 对纹理采样得到的结果就是振幅            sumColumn = getAdd(sumColumn,complexV); // 积分        }        sumRow = getAdd(sumRow,getMultiply(sumColumn, complexU)); //列积分的结果去乘开始那个行上那个欧拉公式算出来的复数, 然后积分    }    float amplitude = sqrt(sumRow.real * sumRow.real + sumRow.imaginary * sumRow.imaginary); // 振幅    float phase = atan2(sumRow.imaginary , sumRow.real); //计算相位, 这里个有坑,atan() 和 atan2() 是两个不同的方法,我第一次用的是atan,后来发现在复数域atan的取值范围不对    rtFourierSpectrum[id.xy] = amplitude/(TextureSize*TextureSize);//clamp(spectrum/TextureSize,0,100); // 这里除以TextureSize是为了好看, 这里放到下面除也可以    rtFourierPhase[id.xy] = phase;    //对两个RT赋值, 你也可以存到一个RT的两个通道里, 放到两个RT上是为了看着方便    rtReal[id.xy] = amplitude * cos(phase);    rtImaginary[id.xy] = amplitude * sin(phase);}

这时候我们就得到了频谱图和相位图,你可以在这里对频谱或者相位进行一些操作。

例如我注释掉那个 //clamp(spectrum/TextureSize,0,0.5); 可以过滤一些高频信息。

或者用各种滤波器在这里处理频谱图, 频谱图代表图片的灰度信息,可以理解为亮度, 而相位代表位置信息。

然后用逆变换还原回去，我们就直接还原回去:

void FourierInverse(uint3 id : SV_DispatchThreadID){    float flag = 2 * 3.14 / TextureSize; // 把e ^ 2π /M 放到这里, 没必要每次循环都算这个    int2 coord = id.xy;//getCoordinate(id.x,id.y);//    Complex sumRow;    sumRow.real = 0;    sumRow.imaginary = 0;    for(int i = 0; i< TextureSize; i++) { // 循环当前行        Complex complexU = getEuler(flag*i*coord.x); //先算当前行首位的的欧拉数 下面列积分后要乘这个数        Complex sumColumn;        sumColumn.real = 0;        sumColumn.imaginary = 0;        for(int j = 0; j < TextureSize; j++) { // 循环每行对应的列            Complex complexV = getEuler(flag*j*coord.y);// 当前像素的欧拉数 - 2 π * v * y / Len            Complex sampleComplex;            //实部虚部、频谱相位是有相互转换关系的，幅度是A，相位是ω的话，实部就是A*cos ω，虚部就是A*sin ω            float amplitude = rtFourierSpectrum[int2(i,j)].r;            float phase = rtFourierPhase[int2(i,j)].r;            sampleComplex.real = amplitude * cos(phase); // 根据公式 a * cos(Phase)  还原实部            sampleComplex.imaginary = amplitude * sin(phase);; // 根据公式 a * sin(Phase)  还原回去虚部            complexV = getMultiply(complexV,sampleComplex); //和该位置的欧拉数相乘            sumColumn = getAdd(sumColumn,complexV); // 积分        }        sumRow = getAdd(sumRow,getMultiply(sumColumn, complexU)); //列积分的结果去乘开始那个行上那个欧拉公式算出来的复数    }    float value = sqrt(sumRow.real * sumRow.real + sumRow.imaginary * sumRow.imaginary); // 求模长 就是求振幅    rtFourierInverse[id.xy] = value;}

最后给一点tips:

• 原始变换公式和逆变换公式其实是完全可以脱离频域和相位的，所以可以直接把复数的实部虚部分开存到两个通道里，让变换流程正常,验证程序没问题
• 频域其实就是复数的模长，相位就是atan(i/r)，所以频域和相位逆变换就是利用三角函数还原复数
• 如果只用复数做变换，那其实是没有意义的，有个频域图和相位图，才有进一步处理合成的空间
• atan() 的取值是  , atan2()的取值是  , 在频域中,我们显然需要的是后者,我开始用的是前者, 卡了好久, 以为算法有问题
• 因为变换中几乎是没办法debug的，因为你没办法知道过程中的值到底是多少(当然不排除利用矩阵硬算)，所以有一组已经验证的复数作为参考可以方便很多, 可以用极简图形先验证频谱正常，在做逆变换，很多图像可以直接看出频谱图应该是什么样，比如以下这个图形: