离散二维傅里叶变换
一常用性质：
       可分离性、周期性和共轭对称性、平移性、旋转性质、卷积与相关定理；
（1）可分离性：
   二维离散傅里叶变换DFT可分离性的基本思想是DFT可分离为两次一维DFT。因此可以用通过计算两次一维的FFT来得到二维快速傅里叶FFT算法。根据快速傅里叶变换的计算要求，需要图像的行数、列数均满足2的n次方，如果不满足，在计算FFT之前先要对图像补零以满足2的n次。
   一个M行N列的二维图像f(x,y)，先按行队列变量y做一次长度为N的一维离散傅里叶变换，再将计算结果按列向对变量x做一次长度为M傅里叶变换就可以得到该图像的傅里叶变换结果，如式所示：
                      

将上式分解开来就是如下的两部分，先得到F(x,v),再由F(x,v)得到F(u,v)：
                        



计算过程如下：


每一行由N个点，对每一行的一维N点序列进行离散傅里叶变换得到F(x,u),再对得到F(x,u)按列向对每一列做M点的离散傅里叶变换，就可以得到二维图像f(x,y)的离散傅里叶变换F(u,v).
同样，做傅里叶逆变换时，先对列向做一维傅里叶逆变换，再对行做一维逆傅里叶变换，如下式所示：


（2）周期性和共轭对称性
由傅里叶变换的基本性质可以知道，离散信号的频谱具有周期性。离散傅里叶变换DFT和它的里变换都以傅里叶变换的点数N为周期的。
对于一维傅里叶变换有：
对于二维傅里叶变换有：
类似有：即从DFT角度来看，反变换得到的图像阵列也是二维循环的。
共轭对称性
对于一维信号有：F(u)=F*(-u),如图所示的一维信号的幅度谱：点数为M的傅里叶变换一个周期为M，关于原点对称。原点即为0频率点，从图中可以看出在0频率的值最大，即信号f(x)的直流分量(均值)，远离原点处的即为高频成份，高频成份的幅值较小，说明信号的大部分能量集中在低频部分。


对于二维信号有：F(u,v)=F*(-u,-v)对于二维图像，其结果如图c所示。左上角(0,0)处为二维图像得0频率点，该点得值对应图像的平均灰度值，图中四个角对应低频成分，中间区域为高频成份，低频区域的幅度值打羽高频区域的幅度值，也同样表示该信号的主要能量集中在低频区域。


根据周期性和共轭对称性，在对图像进行频谱分析处理时只需要关注一个周期就可以了，同时利用图像的傅里叶变换和傅里叶变换的共轭可以直接计算图像的幅度谱，因此使得图像的频谱计算和显示得以简化。
（3)平移性：
傅里叶变换对有如下平移性质：


式子表明，
在频域中原点平移到(u0 ,v0)时，其对应的空间域 f(x,y)要乘上一个正的指数项：
                             

在空域中图像原点平移到(x0,y0)时，其对应的F(u,v)要乘上一个负的指数项：
                              

在数字图像处理中，常常需要将F(u,v)的原点移到N*N频域的中心，以便能清楚地分析傅里叶谱的情况，平移前空域、频域原点均在左上方。要做到这点，只需令上面平移公式中的：u0=v0=N/2；


所以
上式表明：如果需要将图像傅里叶谱的原点从左上角(0,0)移到中心点(N/2,N/2),只要f(x,y)乘上因子进行傅里叶变换即可实现。
平移性还体现了：当空域中f(x,y)产生移动时，在频域中只发生相移，并不影响他的傅里叶变换的幅度，因为：
                                        

反之，当频域中F(u,v)产生移动时，相应f(x,y)在空域中也只发生相移，不产生幅值变化。根据平移性质，为了更清楚查看二维图像的频谱，使直流成分出项在图像中央，在把画面分成四分的基础上，进行如图所示的换位(移位)也是可以的，这样，频域原点就回平移到中心。如下所示：



（4）旋转性质

如果 f(x,y)旋转了一个角度，那么 f(x,y)旋转后的图像的傅立叶变换也旋转了相同的角度。平面直角坐标改写成极坐标形式：

替换则有：
如果f(x,y)被旋转W,则F(u,v)被旋转同一角度。即有傅里叶变换对：
                            


如下所示：



同时，我们可以得出结论，对图像进行旋转变换和傅立叶变换的顺序是可交换的。即先旋转再傅里叶变换或者先傅里叶变换再旋转，得到的结果相同。F{R{f(x,y)}} = R{F{f(x,y)}}。
（5）卷积与相关定理
卷积定理包括空间域卷积和频率域卷积，卷积是空间域滤波和频率域滤波之间的纽带：两个空域信号的卷积等价于其频域信号的 乘积f(x,y)*h(x,y) → F(u,v)H(u,v) 或者 F{f(x,y)*h(x,y)} = F(u,v)H(u,v)；
两个信号频域上的卷积等价于空间域的相乘f(x,y) g(x,y) →F(u,v)*H(u,v)；

该性质的好处是将需要经过翻折、平移、相乘、求和等步骤实现的复杂的卷积运算简化为简单的乘法运算，这也是快速傅里叶变换（FFT）的出现使得该性质得到更广泛应用，同时，该性质对于理解信号的频率域处理方法特别重要，使得信号的空间域处理可以转换到频率域进行处理实现。
根据空间域卷积定理，在空间域对应的是原始信号与滤波器的冲击响应的卷积，卷积定义式为信号翻折平移求和的过程，步骤复杂，运算量大，但如果转换到频率域进行处理，则对在将二者的频谱直接相乘就可以得到滤波结果，然后对滤波结果进行傅里叶逆变换就可以得到滤波后的空间域域图像。如下图所示，对信号进行低通和高通滤波处理的过程和效果。


相关定理：
空域中 f(x,y)与 与 g(x,y) 的相关等价于频域中 F(u,v) 的共轭与 G(u,v)  相乘f(x,y) g(x,y) → F*(u,v)G(u,v)
同时有：f*(x,y)g(x,y) → F(u,v) G(u,v)
相关定理与卷积定理类似，也是把积分求和过程转化为了频域相乘，因此，也使得相关分析的计算简化。
相关的重要应用在于匹配：确定是否有感兴趣的物体区域。f(x,y)是原始图像，g(x,y)作为感兴趣的物体或区域（模板），如果匹配，两个函数的相关值会在 f 中找到相应 g 点的位置上达到最大值。如下图所示。图像 f(x,y) 与模板 g(x,y)，通过计算相关函数，在匹配点处达到最大值，如图中红色圆圈标注的区域。




延拓图像 f(x,y)，延拓图像 g(x,y)，相关函数图像，通过相关图像最大值的水平灰度剖面图。

傅里叶变换的实例与应用
首先我们认识几点有关傅里叶变换的特点：
l 傅里叶变换是从将图像从空间域变换到频率域，具有明确的物理意义。图像的频率是表征图像中灰度变化剧烈程度的指标，是灰度在平面空间上的梯度，在噪声点和图像边缘处的频率为高频。
l 在频率域中，将信号表示为一系列正弦信号或者复指数函数的叠加，正弦信号的频率、幅值和相位可以描述正弦信号中的所有信息，由此可以得到信号的幅度谱和相位谱。在图像领域就是将图像灰度作为正弦变量。
l 傅里叶变换全局性的，是一个积分求和的过程，对时间、地点位置无法进行准确定义，也就是说傅里叶变换得到的频谱图中的点无法与空间域中的某个空间位置对应，因此，从傅里叶变换图中并不能直接对应某个位置的特点。
l 傅里叶变换是一系列不同频率三角函数的和，每个频率分量的系数不同，这些系数代表了各频率成分的强弱或者所占比重，通过分析这些系数就可以分析图像的特性。

1、介绍。
 
        DFT:(Discrete Fourier Transform)离散傅里叶变换是傅里叶变换在时域和频域上都呈离散的形式，将信号的时域采样变换为其DTFT的频域采样。在形式上，变换两端（时域和频域上）的序列是有限长的，而实际上这两组序列都应当被认为是离散周期信号的主值序列。即使对有限长的离散信号作DFT，也应当将其看作其周期延拓的变换。
 
1）、欧拉公式：
 
，其中i是虚数，即i的平方为-1。
 
 
 
2）、二维离散傅里叶变换DFT公式：
 
N是二维数组的行数，M是二维数组的列数。u和v是转换后二维数组的位置，F(u,v)是转换后数组中相应位置的值。x和y是原二维数组的位置，f(x,y)是原数组中相应的值。
 

快速傅里叶变换 (fast Fourier transform), 即利用计算机计算离散傅里叶变换（DFT)的高效、快速计算方法的统称，简称FFT。快速傅里叶变换是1965年由J.W.库利和T.W.图基提出的。采用这种算法能使计算机计算离散傅里叶变换所需要的乘法次数大为减少，特别是被变换的抽样点数N越多，FFT算法计算量的节省就越显著。
 具体的贴出学习资源，个人认为讲解的非常好。
 https://blog.csdn.net/ggn_2015/article/details/68922404
 click
 click
 知识点以及证明综合都有利于理解。
 
上一道例题Gym100783C
 题意：给定n个数的数组a，m个数的数组b，问b数组中的数有多少个可以由a数组中的1或者2个数值和组成，可以重复在a数组中选取。
 
可以看成，两个多项式的相乘，指数级别是相加的，系数分存在或者不存在，即0或者1。即a数组*a数组，看最后中的指数部分是否在b中出现。
 
自己非常菜，改的板子。
 
#include<cmath>
#include<iostream>
#include<cstdio>
#include<cstring>
#include<string>
#include<cstdlib>
#include<istream>
#include<vector>
#include<stack>
#include<set>
#include<map>
#include<algorithm>
#include<queue>
#define inf 0x3f3f3f3f
#define llinf 0x3f3f3f3f3f3f3f3f
#define MAX_len 210*4
using namespace std;
typedef long long ll;
const int MAXN = 262144 * 2+ 10;
const double PI = std::acos(-1.0);

struct Complex {
    double r, i;

    Complex(double r = 0, double i = 0) : r(r), i(i) {}

    Complex conj() const { return Complex(r, -i); }

    Complex operator+(const Complex &rhs) const { return Complex(r + rhs.r, i + rhs.i); }
    Complex operator-(const Complex &rhs) const { return Complex(r - rhs.r, i - rhs.i); }
    Complex operator*(const Complex &rhs) const { return Complex(r * rhs.r - i * rhs.i, r * rhs.i + i * rhs.r); }
    Complex operator/(double rhs) const { return Complex(r / rhs, i / rhs); }
};

class FFT {
private:
    static const int N = 262144*2;

    Complex omega[N + 1], omegaInv[N + 1];

    void init() {
        for (int i = 0; i < N; i++) {
            omega[i] = Complex(std::cos(2 * PI / N * i), std::sin(2 * PI / N * i));
            omegaInv[i] = omega[i].conj();
        }
    }

    void reverse(Complex *a, int n) {
        for (int i = 0, j = 0; i < n; i++) {
            if (i < j) std::swap(a[i], a[j]);
            for (int l = n >> 1; (j ^= l) < l; l >>= 1) {}
        }
    }

    void transform(Complex *a, int n, Complex *omega) {
        reverse(a, n);

        for (int l = 2; l <= n; l <<= 1) {
            int hl = l >> 1;
            for (Complex *x = a; x != a + n; x += l) {
                for (int i = 0; i < hl; i++) {
                    Complex t = omega[N / l * i] * x[i + hl];
                    x[i + hl] = x[i] - t;
                    x[i] = x[i] + t;
                }
            }
        }
    }

public:
    FFT() { init(); }

    int extend(int n) {
        int res = 1;
        while (res < n) res <<= 1;
        return res;
    }

    void dft(Complex *a, int n) {
        transform(a, n, omega);
    }

    void idft(Complex *a, int n) {
        transform(a, n, omegaInv);
        for (int i = 0; i < n; i++) a[i] = a[i] / n;
    }
} fft;
int aaa[MAXN];
map<int,int>hh1;
int main() {
    int n, m;
    scanf("%d", &n);

    static Complex a[MAXN], b[MAXN];
    int MAXXX=0;
    for(int i=0;i<n;i++)
    {
        scanf("%d",&aaa[i]);
        hh1[aaa[i]]=1;
        MAXXX=max(MAXXX,aaa[i]);
    }
    for (int i = 0, x; i <=MAXXX; i++) {
        if(hh1[i])
        a[i] = Complex(1, 0);
        else
             a[i] = Complex(0, 0);
    }

    for (int i = 0; i <=MAXXX; i++) {
        if(hh1[i])
        a[i] = a[i] + Complex(0,1);
        else
                    a[i] = a[i] + Complex(0,0);
    }

    int len =MAXXX + MAXXX + 1;
    int N = fft.extend(len);

    fft.dft(a, N);
    for (int i = 1; i < N; i++) {
        double x1 = a[i].r, y1 = a[i].i;
        double x2 = a[N - i].r, y2 = a[N - i].i;
        Complex t1((x1 + x2) * 0.5, (y1 - y2) * 0.5);
        Complex t2((y1 + y2) * 0.5, (x2 - x1) * 0.5);
        b[i] = t1 * t2;
    }
    b[0] = a[0].r * a[0].i;
    fft.idft(b, N);
    int ans=0;
    for (int i = 0; i < len; i++)
    {
        int temp=(round(b[i].r));
        if(temp)
        {
            hh1[i]=1;
        }
    }
    scanf("%d",&m);
    for(int i=0;i<m;i++)
    {
        int x;
        scanf("%d",&x);
        if(hh1[x])
            ans++;
    }
     printf("%d",ans);

    return 0;
}
/*
3
1 3 5
6
2 4 5 7 8 9
*/

首先安利一个网站，在线做傅里叶变换，不用等MATLAB漫长的启动了
 https://sci2fig.herokuapp.com/fourier
 文章中部分图片来自
 https://blog.csdn.net/ViatorSun/article/details/82387854
 理论部分参考了：
 https://homepages.inf.ed.ac.uk/rbf/HIPR2/fourier.htm
 https://plus.maths.org/content/fourier-transforms-images
 


 
零基础也能理解傅里叶变换后的频谱图
 
频谱图
类比
  
图像的傅里叶变换原理
位置信息
图像高频低频含义
  
对图像的操作所造成的影响
具体公式和理论应用
常见的傅里叶变换结果
练习题
解答
 

 
频谱图
 
声音的频谱图很好理解，尖峰代表着该频率有着更多的分量
 
 关于数字信号处理的更多内容，这里安利我的学习笔记【持续更新中】：【从零开始学信号与系统】
 
但是图像的傅里叶变换抽象的多，下面讲讲述如何看懂傅里叶图。
 
类比
 
将图像视为变化的函数，不过不是随时间变化，而是跨图像的二维空间变化。
 在灰度数字图像中，每个像素的值在0到255之间，代表该像素的暗度。因此，该像素的暗度或强度是水平坐标和垂直坐标的函数，给出了该像素的位置。您可以将图像视为起伏的景观，而高度由像素值确定。
 

 图像也可以表示为正弦波的总和，但是这次，它们不是一维波，而是二维变化的波，就像纸张上的波纹。
 
二维正弦波写为
 
$z=sin（hx+ky）$
 
其中x和y给出“图纸”上各点的坐标，z是该点处波浪的高度或强度，a给出波浪的振幅（最大高度），h和k给出数字波分别在x和y方向重复的次数（它们是 x和y频率）。
 
当k = 0时，正弦波仅沿x轴波动。当h = 0时，它仅沿y轴波动。但是，如果k 和h都不为零，则正弦波在片材上对角移动，并且波的传播方向（垂直于波阵面）与斜率h / k成一定角度 。
 
图像的傅里叶变换原理
 
图像的傅立叶变换将图像（起伏的地形）分解为正弦波的总和。就像声波一样，对频率绘制傅立叶变换。但是与那种情况不同，频率空间具有二维，对于x和y维 中的波的频率 h和k。因此，它不是以一系列尖峰的形式绘制的，而是以与原始图像（大约）相同的像素尺寸绘制的图像。
 
位置信息
 
傅立叶变换中的每个像素都有一个坐标（h，k）表示在傅立叶变换中具有x频率h和 y频率k的正弦波的贡献。
 
中心点表示（0,0）波–没有波纹的平面–其强度（灰度颜色的亮度）是图像中像素的平均值。
中心左侧和右侧的点代表沿x轴变化的正弦波，即k = 0）。这些点的亮度表示傅立叶变换中具有该频率的正弦波的强度（强度是正弦波的振幅的平方）。
在中心点上下垂直的那些代表那些在y中变化但在x中保持恒定（即h = 0）的正弦波。傅立叶变换中的其他点表示对角波的贡献。
 
 例如，考虑上面的图像，在左侧。这是二维波sin（x）被视为灰度图像。它的旁边是此灰度图像的傅立叶变换。
它具有与原始像素相同的尺寸（像素），并且全黑
在中心处有一些明亮的像素
放大傅立叶变换的中心（您可以在上面在右侧看到），您会看到恰好有三个不是黑色的像素。 
  
一个是明亮的中心点，坐标为（0,0），代表（0,0）波对图像的贡献。
两侧的明亮像素具有坐标（1,0）和反射（-1,0），表示（1,0）波（我们原始图像中的正弦波）的贡献。
傅立叶变换中的所有其余像素都是黑色的，因为仅使用原始（1,0）波精确描述了原始图像。
 
 
 
图像高频低频含义
 
图像的频率是表征图像中灰度变化剧烈程度的指标，是灰度在平面空间上的梯度。
 
简短概括：
 
图像的高频，意味着灰度变化剧烈
图像的低频，意味着灰度变化平坦
 
不同频率信息在图像结构中有不同的作用：
 
图像的主要成分是低频信息，它形成了图像的基本灰度等级，对图像结构的决定作用较小；
中频信息决定了图像的基本结构，形成了图像的主要边缘结构；
高频信息形成了图像的边缘和细节，是在中频信息上对图像内容的进一步强化。
 
图像进行二维傅立叶变换得到频谱图，就是图像梯度的分布图,当然频谱图上的各点与图像上各点并不存在一一对应的关系，即使在不移频的情况下也是没有。
 傅立叶频谱图上我们看到的明暗不一的亮点，实际是上图像上某一点与邻域点差异的强弱，即梯度的大小，也即该点的频率的大小（可以这么理解，图像中的低频部分指低梯度的点，高频部分相反）。
 
对图像的操作所造成的影响
 
在二维傅里叶变换中，空间域中横向的周期变化会反应在频谱图中的Y轴上，而空间域中纵向的周期变化会反应在频谱图中的X轴上。空间域中东南方向的周期变化会反应在频谱图中的东北方向，反之亦然。说明见下图。
 
 
 
具体公式和理论应用
 
参考：https://homepages.inf.ed.ac.uk/rbf/HIPR2/fourier.htm
 由于我们只关注数字图像，因此我们将讨论仅限于离散傅立叶变换（DFT）。
 
DFT是采样的傅立叶变换，因此不包含构成图像的所有频率，而仅包含足够大以完全描述空间域图像的一组采样。频率的数量对应于空间域图像中的像素数量，即，空间域和傅立叶域中的图像具有相同的大小。
 
对于大小为N×N的正方形图像，二维DFT由下式给出：
 
 其中f（a，b）是空间域中的图像，指数项是与傅立叶空间中每个点F（k，l）对应的基函数。该方程式可以解释为：通过将空间图像乘以相应的基函数并将结果相加得出每个点F（k，l）的值。
 
基本函数是频率增加的正弦和余弦波，即 F（0,0）表示图像的DC分量，该分量对应于平均亮度，F（N-1，N-1）表示最高频率。
 
以类似的方式，可以将傅立叶图像重新变换到空间域。傅里叶逆变换由下式给出：
 【这里不太会翻译，看图吧】
 
 原网站还有关于低通滤波等操作的详细知识，吃透了再补文。
 
常见的傅里叶变换结果
 
不同频率的正弦波：
 
 MIT的摄影师：
 
 
练习题
 
判断上面两张傅里叶变换频谱图和下面两张原始图像如何对应？
 
 
解答
 
右下角的图片灰度变化更剧烈，由于图像的频率是灰度的梯度，且傅里叶变换后高频特征位于边缘，因此对应着亮点分布在更边缘地区的左上图：