快速傅里叶变化：主要有几个难点

1.如果使用DIT首先要倒位序，网上的方法值得参考

思路是和二进制进位的思路相同只是是反着进位

2.要自己封装复数类型

3.难点是公式里的k不好求，因为k不是连续取得的，我用了一个数j来存储，这样可以把k计算出来

4.wn^r中r的取值，直接移位即可，不需要转化为二进制

5.输入不是2的倍数时自动补零

其余我都写在注释里了，感觉还是有很多地方要改进
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include "复数类型.c"
//可以用M_E表示自然对数，M_PI表示圆周率

/*首先要导位序，也可以只要一个存储空间A[i],只需要改变A[i]中i
重点1.只改变当i小于j的时候的值，而i大于等于j时因为已经换过，不能再换回来了
重点2.思路是与二进制算法思路，小于N/2时即，最高位为0，此时的结果需要加上N/2
当大于N/2时，最高位为1，此时要向低位进位，意味着要先-N/2代表进位
再判断是否大于N/2，大于再向次低位进位，否则直接加N/4*/

void daoxu(double *x,int n)
{
    int i=0,j=0,k=0;//i表示原始值，而j表示倒位序后的值
    double temp;
    for(;i < n;i++)
    {
        if(i<j)
        {
            temp = x[j];
            x[j] = x[i];
            x[i] = temp;
        }

        k = n/2;
        while(j >= k && k != 0)
        {
            j = j - k;
            k = k/2;
        }
        j = j + k;     
    } 
    /*for(i=0;i<n;i++)
    {
        printf("x[%d]=%f\n",i,x[i]);
        //printf("result: new_k[%d]=%d\n",i,new_k[i]);
    }*/
}

int process(int* ori_len)//对输入的数据进行处理，还要对数据进行补零，问题在于如何分配内存，使用链表？
{
    int temp = 1,i;
    for(i=0;;i++)
    {
        if(temp >= *ori_len)
            break;
        temp =temp * 2;
    }
    *ori_len = temp;
    return i;//返回幂次方
}

complex w_n(int k,int m,int L,int length)
//如果这个值是复数意味着要包装复数类型，并且和这个数做四则运算的其他数也必须是复数
//需要新建复数这个数据结构，但是其实最后计算的是幅值
{
    int r=k<<(L-m);//不应该取这个k，而应该取倒位序的k
    r=r&(length-1);//去除高位
    //printf("得到的r为=%d\n",r);
    complex result;
    result.real=cos(2*M_PI/length*r);//用欧拉公式展开
    //printf("result.real=cos(2*pi/%d*%d)=%f\n",length,r,result.real);
    result.imag=-sin(2*M_PI/length*r);
    //printf("result.imag=-sin(2*pi/%d*%d)=%f\n",length,r,result.imag);
    return result;
}

void butterfly_operation(double *x,int powe,int length)
{
    int k_num = length/2,k_cal,k,l;
    complex r;
    complex *back_up=(complex*)calloc(length,sizeof(complex));//需要备份值来保留处理前的数据,定义为虚数类型
    complex *y=(complex*)calloc(length,sizeof(complex));
    
    for(int m=1;m<=powe;m++)//pow表示蝶形的次数
    {
        //printf("第%d次蝶形运算\n",m);
        //printf("要做%d次蝶形运算\n",k_num);
        k = 0;//运算之前将k清零,这里的k是实际是倒位序的
        l = 0;
        for(int i=0;i<length;i++)//运算完一次蝶形更新备份值
        {
            if(m==1)
                back_up[i] = cpx_make(x[i],0);
            else
                back_up[i] = y[i];
            //printf("back_up[i]=%f+%fi\n",back_up[i].real,back_up[i].imag);
        }
        k_cal=pow(2,m-1);//k_cal实际代表的是2^(m-1)的值
        //printf("计算2^(m-1)的值=%d\n",k_cal);
        for(int j=0;j<k_num;j++)//在做第m次蝶形的时候，要做k_num次蝶形运算
        {
            //printf("k=%d\n",k);
            r = w_n(k,m,powe,length);
            //printf("r=%f+%fi\n",r.real,r.imag);
            y[k] =cpx_add(back_up[k],cpx_mul(back_up[k+k_cal],r));//x[k]取的是模值
            //printf("y[%d]=%f+%fi\n",k,y[k].real,y[k].imag);
            y[k+k_cal] = cpx_sub(back_up[k],cpx_mul(back_up[k+k_cal],r)); 
            //printf("cpx_mul(back_up[k+k_cal],r)=%f+%fi\n",cpx_mul(back_up[k+k_cal],r).real,cpx_mul(back_up[k+k_cal],r).imag);
            //printf("y[%d]=%f+%fi\n",k+k_cal,y[k+k_cal].real,y[k_cal+k].imag);
            if(k_cal>=2&&k<l+k_cal-1)
            {
                k+=1;               
            }
            else
            {
                k=k+k_cal+1;
                l=k;
            }      
        }
        //printf("-------------------------------------------------------------------\n");
    }
    for(int i=0;i<length;i++)
    {
        x[i]=sqrt(y[i].real*y[i].real+y[i].imag*y[i].imag);
        //printf("y[i].real=%f,y[i].imag=%f\n",y[i].real,y[i].imag);
    }
    free(back_up);
    free(y);
}

int main()
{
    int length,power,old_length;//以后要更正
    double t;
    printf("请输入采样的长度length=");
    scanf("%d",&length);//输入采样的长度*/
    double *x,*y;
    y=(double*)calloc(length,sizeof(double));//可以动态分配内存
    for(int i=0;i<length;i++)
    {
        //y[i] = i+1;//在输入之前就会分配好内存空间因此不能临时输入
        t=i*0.001;
        y[i] = 1.2+2.7*cos(2*M_PI*33*t)+5*cos(2*M_PI*200*t+M_PI/2);
        //y[i]=0.6*sin(2*M_PI*500*i)+0.6*sin(2*M_PI*50*i);
        //y[i] = sin(M_PI*i/8);
        //printf("%f\n",y[i]);
    }
    old_length = length;
    power = process(&length);//改变了length    
    x=(double*)calloc(length,sizeof(double));   
    for(int i=0;i<old_length;i++)
    {
        x[i] = y[i];
        printf("x[i]=%f\n",x[i]);
    }
    for(int i=old_length;i<length;i++)
    {
        x[i]=0;
        printf("x[i]=%f\n",x[i]);
    }
    free(y);
    daoxu(x,length);
    butterfly_operation(x,power,length);
    for(int i=0;i<length;i++)
        printf("%d\t%f\n",i,x[i]);//在输入之前就会分配好内存空间因此不能临时输入
    free(x);
    return 0;
}

我自己封装的复数头文件
#include <stdio.h>
#define Real(c) (c).real
#define Imag(c) (c).imag

typedef struct
{
    double real;
    double imag;
    
}complex;

complex cpx_make(double real, double imag)
{
    complex ret;
    ret.real = real;
    ret.imag = imag;
    return ret;
}

complex cpx_add(complex a,complex b)
{
    return cpx_make(Real(a) + Real(b), Imag(a) + Imag(b));
}

complex cpx_sub(complex a,complex b)
{
    return cpx_make(Real(a) - Real(b), Imag(a) - Imag(b));
}

complex cpx_mul(complex a,complex b)
{
    return cpx_make(Real(a) * Real(b)-Imag(a) * Imag(b),Real(a) * Imag(b) + Imag(a) * Real(b));
}

complex cpx_div(complex a,complex b)
{
    Real(b) = 1/Real(b);
    Imag(b) = 1/Imag(b);
    return cpx_mul(a,b);
}

最早的电话使用的模拟信号原理图

（1）声音通过金属振动膜感应声波来影响磁场和电流，并将这种带有金属振动膜振动的“信息”的电流传递给另一端
（2）另一端则进行反向工作，把不断变化的电流转化为电线圈中磁场的变化，使得金属振动膜发出同样的振动。
局限性就是当距离比较远 的时候就很难还原（电流的衰减大）
傅里叶变换和反傅里叶变换就解决了这个问题--使得数字信号通信代替模拟信号通信称为可能
（1）任意一个是与信号f(t)经过傅里叶变换，都可以转化为多个余弦波叠加的形式
一旦这个环节被我们掌控的话，语音信号就可以通过调制，由高频载波用极高的速率发送电信号和光信号。
而在接收端，用滤波器过滤余弦波信号后，使用傅里叶反变换进行是与信号的还原

对傅里叶的理解资源：
深入浅出的讲解傅里叶变换（真正的通俗易懂） - CSDN博客
https://blog.csdn.net/l494926429/article/details/51818012
实验的来源：
实验六傅里叶变换及其反变换_百度文库
https://wenku.baidu.com/view/1ba79654178884868762caaedd3383c4bb4cb4bd.html

代码：
syms t v w x phase im re ;          %定义符号变量  
a = input('请输入a=')
f = exp(-a*abs(t));  %f(t) = exp(-2*t)*u(t)
Fw = fourier(f);                     %求傅里叶变换
subplot(311);
ezplot(f);                          %绘制f（t）的时域波形
axis([-1 2.5 0 1.1]);
subplot(312);
ezplot(abs(Fw));                    %绘制幅度谱
im = imag(Fw);                      %计算F(w)的虚部
re = real(Fw);                      %计算F(w)的实部
phase = atan(im/re);                %计算相位谱
subplot(313);
ezplot(phase);                      %绘制相位谱

syms t v w x phase im re ;          %定义符号变量  
F(w) = 1/(1+w^2);                  %待求解变换的公式
f = fourier(Fw,t);                 
subplot(311);
ezplot(f);                          %绘制f(t)的时域波形
axis([-1 2.5 0 1.1]);
subplot(312);
ezplot(abs(Fw));                    %绘制幅度谱
im = imag(Fw);                      %计算F(w)的虚部
re = real(Fw);                      %计算F(w)的实部
phase = atan(im/re);                %计算相位谱
subplot(313);
ezplot(phase);                      %绘制相位谱

快速傅里叶变化（FFT）是离散傅里叶变化（DFT）的快速算法。一个长度为N的有限长序列的DFT为

式中的 为 。离散傅里叶反变换为

FFT算法的基本思想是利用 的对称性和周期性，即 =- 和 （r为任意整数），通过将长序列分解为短序列，再由短序列的DFT组合来实现整个长序列的DFT，从而减少DFT的复数乘法和复数加法次数，提高DFT的运算速度和效率。将N点的长序列分解为短序列可以有多种方式，常用的有基2时间抽取FFT算法，基2频率抽取FFT算法，基4抽取FFT算法和分裂基FFT算法。基2抽取算法简单有效，基4抽取算法和分裂基抽取算法与基2抽取算法原理类似，但基4算法和分裂基算法效率更高，运算量更少。本文主要介绍基2抽取FFT算法。
   1、基2时间抽取FFT算法：
基2时间抽取FFT算法是将N点序列x[n]按序号奇偶分解得到两个 点子序列 和 ,通过计算这两个子序列的DFT 和 实现 的DFT 。算法要点为N点序列x[n]对时域下标按奇数和偶数分解，即
（1）

对频域X[k]按前一半和后一半分解，即 。频域的前一半可由下式求得

由DFT的周期性和 的对称性有


-=
故频域的后一半可由下式求得

将前一半和后一半组合起来，即得原来的长序列x[n]的DFT 。
可对 和 继续按式(1)分解下去，直到每个子序列的长度为1为止。任何一个N= 点的DFT，可由通过M次分解，得到长度为1的子序列。对于N不是2的整数次幂的情况，可在序列后面添0补齐到2的整数次幂。每级分解运算可用一个蝶形运算流程图来表示，如图

下图是对一个长度为8的序列进行基2时间抽取FFT运算的流程图

直接计算N点序列的DFT需要N2次复数乘法，需要N(N-1)次复数加法。而用基2时间抽取FFT算法，则乘法次数为 ，加法次数 ，运算次数大大减少。
2、基2频率抽取FFT算法
基2频率抽取FFT算法则是对频域下标按奇数和偶数分解，对时域下标按前一半和后一半分解。 的DFT 可表示为

将 分解成 和 两个长度为 的子序列，如下：
（2）
则 的DFT 可按频率下标分奇数和偶数表示为


可继续将 和 按式(2)分解下去，直到每个子序列长度为1。任何一个N= 点的DFT，可由通过M次分解，得到长度为1的子序列。对于N不是2的整数次幂的情况，可在序列后面添0补齐到2的整数次幂。每级分解运算也可用一个蝶形运算流程图来表示，如图

下图是对一个长度为8的序列进行基2频率抽取FFT运算的流程图：

3、编程实现
Matlab有一个内置的函数fft()可以实现快速傅里叶变换，由于它是built-in函数，所以无法查看源码来判断它是应用的那种分解方式来实现的快速傅里叶变换。下面将用Matlab语言来编写实现类似fft()的快速傅里叶变换函数myfft()。本函数采用的FFT算法为基2时间抽取FFT算法，函数定义源码如下：
%----------------------------------
%myfft(A,M)有两个形参，A为要进行快速傅里叶变换的序列，M为序列长度对2的对数，
%即如果序列A[n]的长度为8，则M=3
function [A] = myfft(A,M)
N=2^M;                   %序列的长度为N
LH=N/2;                  %序列长度的一半
N1=N-2;
J=LH
for I=1:1:N1             %将输入序列A[n]按运算流程图第一列的顺序排好
    if I<J
        temp=A(I+1);
        A(I+1)=A(J+1);
        A(J+1)=temp;
    end
    K=LH;
    while J>=K
        J=J-K;
        K=K/2;
    end
    J=J+K;
end
for L=1:1:M              %L表示序列分解层次，此次分解有M-L个子序列
    B=2^(L-1);
    for S=0:B-1          %执行第S个蝶形运算，每个子序列共有B个蝶形运算
        p=S*2^(M-L);      %旋转因子的上标，即此时旋转因子为 
        for k=S:2^L:N-1   %执行各个子序列中的第S个蝶形运算，有M-L个子序列
            temp=A(k+1)+A(k+B+1)*exp(-i*2*pi*p/N);
           A(k+B+1)=A(k+1)-A(k+B+1)*exp(-i*2*pi*p/N);
            A(k+1)=temp;  %蝶形运算
        end
    end
end
end
    python语法简洁清晰，具有强大和丰富的三方库，能把其他语言(如C/C++,Matlab，java等)制作的各种模块轻松联结在一起，开发效率高，很适合快速开发实验原型以验证算法正确性和有效性。下面是基2时间抽取FFT算法的一个python版本实现。
#-------------------------------------
import numpy as np
def  FFT(x): 
    x =np.asarray(x, dtype=float)
    N =x.shape[0]
    ifN%2>0:                                                      #输入序列的长度必须为2的幂
        raiseValueError("size of x must be a power of 2")
    elif N<= 2:                                                   # 子序列长度小于等于2时，直接计算
         x =np.asarray(x, dtype=float)
            N = x.shape[0]
            n = np.arange(N)
            k = n.reshape((N, 1))
            M = np.exp(-2j * np.pi * k * n / N)
            eturn np.dot(M, x)
    else:
        X_even= FFT(x[::2])                               #偶序列继续分解
        X_odd= FFT(x[1::2])                                              #奇序列继续分解
        factor= np.exp(-2j * np.pi * np.arange(N/2) / N)        #旋转因子
        returnnp.concatenate([X_even + factor* X_odd,         #返回序列频域前半部分
                            X_even - factor*X_odd])     #返回序列频域后半部分
4、函数测试和实验
快速傅里叶变换可以得到信号的频谱。通过对一个给定的信号进行快速傅里叶变换，观察它的频谱，并和Matlab内置函数fft()变换得到的结果进行对比，来对myfft()的功能进行测试。测试可用如下脚本来实现：
%----------------------------
%测试脚本，分别用myfft()和fft()对一个信号进行快速傅里叶变换，该信号由一个
%50Hz和120Hz的正弦信号叠加，频谱会在这50Hz和120Hz出现峰值。
Fs = 1000;                    % 采样频率
T = 1/Fs;                     % 采样周期
L = 5000;                     % 信号长度
t = (0:L-1)*T;               % 时间序号数组
% x(t)=0.7sin(2 *50t)+sin(2 *150t),为50Hz和150Hz正弦信号的叠加
x = 0.7*sin(2*pi*50*t) + sin(2*pi*150*t); 
y = x + 2*randn(size(t));     % 给信号加上噪声
%画出x(t)和y(t)的图像
subplot(2,1,1)
plot(Fs*t(1:50),x(1:50))
title('Sinusoids Signal')
xlabel('time (milliseconds)')
subplot(2,1,2)
plot(Fs*t(1:50),y(1:50))
title('Signal Corrupted with Zero-Mean Random Noise')
xlabel('time (milliseconds)')
%对加了噪声的y(t)进行快速傅里叶变换
NFFT = 2^nextpow2(L);       %延长信号长度为最小的2的幂
y=[y zeros(1,2^nextpow2(L)-L)]; %延长信号长度，后面用0补齐
Y = myfft(y,log2(NFFT))/L; % 用自己实现的快速傅里叶变化myfft处理信号得到Y（f）
Y2 = fft(y,NFFT)/L;        %用Matlab官方实现的fft处理信号得到Y2(f）
f = Fs/2*linspace(0,1,NFFT/2+1);%频域横坐标
%对信号y(t)进行快速傅里叶变换，并画出的频谱
figure
subplot(2,1,1)
plot(f,2*abs(Y(1:NFFT/2+1)))       %                              

title('myfft()处理的y(t)频谱')% y(t)')
xlabel('Frequency (Hz)')
ylabel('|Y(f)|')
subplot(2,1,2)
plot(f,2*abs(Y2(1:NFFT/2+1)))       %                               
title('fft()处理的y(t)频谱')
xlabel('Frequency (Hz)')
ylabel('|Y2(f)|')
%对原信号x(t)进行快速傅里叶变换，并画出频谱
x=[x zeros(1,2^nextpow2(L)-L)]; %延长信号长度，后面用0补齐
X = myfft(x,log2(NFFT))/L;     % 用自己实现的快速傅里叶变化myfft处理信号得到X（f）
X2 = fft(x,NFFT)/L;             %用Matlab官方实现的fft处理信号得到X2(f）
figure
subplot(2,1,1)
plot(f,2*abs(X(1:NFFT/2+1)))       %                              

title('myfft()处理的x(t)频谱')%
xlabel('Frequency (Hz)')
ylabel('|X(f)|')
subplot(2,1,2)
plot(f,2*abs(X2(1:NFFT/2+1)))       %                              

title('fft()处理的x(t)频谱')
xlabel('Frequency (Hz)')
ylabel('|X2(f)|')
运行脚本，得到x(t)和y(t)的频谱，如图。可见，myfft()和Matlab内置函数fft()的处理结果一致，没有区别。频谱图在50Hz和150Hz出现峰值，幅值分别为0.63和1.0左右。

  
Python版的测试脚本如下
a=np.array([1,0,1,0,1,0,1,0])           #定义序列a[n]=[1,0,1,0,1,0,1,0]
b=[1,1,1,1,0,0,0,0]                         #定义序列b[n]=[1,1,1,1,0,0,0,0]
A=FFT(b)                                       #对a[n]进行快速傅里叶变换
B=FFT(a)                                      #对b[n]进行快速傅里叶变换
printA,”\n”,B                                  #输出两个序列的傅里叶变换的结果A[M]和B[m]
python版的傅里叶变换结果和Matlab自带的fft()变换结果对比如下

Python版的运行结果

Matlab版的运行结果
注意python版本的运行结果中，有”e-16”的结果，这是因为计算机浮点数的运算误差导致的，实际上可认为是0。可见，python实现的快速傅里叶变换函数的结果和Matlab的fft()变换结果一致，python实现的快速傅里叶变换结果正确无误。
5、快速傅里叶变换的应用
快速傅里叶变换可以将有限长或者周期离散时间信号从时域表示转换到频域表示。快速傅里叶变换除了可以得到一个信号的频谱外，还可以基于得到的频谱对信号进行一些频域上的处理，并通过快速傅里叶反变换得到经过处理后的信号。比如计算离散时间信号的线性卷积，滤波，信号压缩，降噪，谱估计等。FFT除了用于离散时间信号处理，也可以用于逼近连续时间信号的频谱。
下面脚本展示了利用快速傅里叶变换和反变换进行快速卷积计算
subplot(3,1,1);
n=0:16;
x=0.8.^n;      %x[n]=0.8^n的离散信号
stem(n,x);
xlabel('n');ylabel('x[n]');
subplot(3,1,2);
n=0:15; 
y=[ones(1,10) zeros(1,6)]; 
stem(n,y);
xlabel('n');ylabel('y[n]')
subplot(3,1,3);
L=26;n=0:L-1;
X=fft(x,L);
Y=fft(y,L);
Z=X.*Y;      %在频域两信号的频谱相乘
z=ifft(Z,L);  %快速傅里叶反变换得到两信号的卷积
stem(n,z); 
xlabel('n');
ylabel('z[n]') 
脚本运行结果如图

在工程领域，FFT可以用于语音处理，工程传感器信号处理，特征提取等。
二维FFT和IFFT可应用到图像处理领域，比如图像降噪，模糊化，清晰化，轮廓化，压缩，还原等。快速傅里叶变换及其反变换是图像处理领域的重要处理手段，是进行其他各项处理的基础，应用场合很广。下面脚本展示了图像的二维快速傅里叶变换及其反变换，及基于快速傅里叶变换得到的频谱图，相位图等一系列频域信息。
C=imread('2.jpg');%载入图片
A=rgb2gray(C);
B=fftshift(fft2(A));% 进行傅立叶变换
E=fft2(A);
D=ifft2(E);
subplot(331)
imshow(C);
title('原始图像');
subplot(332)
imshow(A);
title('原始灰度图像'); 
subplot(339)
imshow(abs(B),[]); 
title('原始频谱图'); 
subplot(334)
imshow(log(abs(B)),[]);
title('取对数后的频谱图'); 
subplot(335)
imshow(angle(B),[]); 
title('相位图'); 
subplot(336)
imshow(real(B),[]); 
title('实部图'); 
subplot(337)
imshow(imag(B),[]); 
title('虚部图'); 
subplot(338)
imshow(abs(E),[]); 
title('fft2频谱图');
subplot(333)
imshow(D,[]); 
title('经快速傅里叶逆变换得到的灰度图');
下面为脚本运行的结果，可见经快速傅里叶反变换后，可以从图像的频域表示得到时域表示，即重新得到一张图像。由于此脚本没有对图像做其他处理，所以快速傅里叶反变换得到的图形与原来的灰度图像一致。

6、总结
快速傅里叶变换大大简化了运算量，并且算法具有迭代和分治的特点，特别适合计算机实现，使得用计算机对离散时间信号进行快速处理成为现实。配合采样定理和信号采样及重构等方法，使得快速傅里叶变换可以进一步推广到连续时间信号处理，大大提高了信号处理的效率和速度。在语音处理领域，傅里叶变换可以用于语音识别，音频分离，降噪等应用。在图像应用领域，快速傅立叶变换的物理意义是将图像的灰度分布函数变换为图像的频率分布函数，傅立叶逆变换是将图像的频率分布函数变换为灰度分布函数。而灰度分布函数即为灰度图的表达。快速傅里叶变换使得对图像可以用频域处理的方法来处理，对基于计算机的图像处理应用的可行性和实时性具有重大意义，是数字信号处理的基石。

目的：读取图像 A(lena.tiff)和B(rice.tif)，显示这两幅图像，对图像作傅立叶变换，显示图像的傅里叶幅度谱和相位谱。做傅立叶逆变换，显示重建图像。

图像的频率是表征图像中灰度变化剧烈程度的指标，是灰度在平面空间上的梯度

对图像而言，图像的边缘部分是突变部分，变化较快，因此反应在频域上是高频分量；图像的噪声大部分情况下是高频部分；图像平缓变化部分则为低频分量。也就是说，傅立叶变换提供另外一个角度来观察图像，可以将图像从灰度分布转化到频率分布上来观察图像的特征。

imshow()函数： 

imshwo()函数用于接收一个像素矩阵，显示该图像，其显示的参数有两种类型 

  unit8;像素在矩阵处理范围为0-255 

  double：若值大于1，转化为1，若小于1，转化为0

图像进行二维傅立叶变换得到频谱图，就是图像梯度的分布图,当然频谱图上的各点与图像上各点并不存在一一对应的关系，即使在不移频的情况下也是没有。傅立叶频谱图上我们看到的明暗不一的亮点，实际是图像上某一点与邻域点差异的强弱，即梯度的大小，也即该点的频率的大小（可以这么理解，图像中的低频部分指低梯度的点，高频部分相反）。 

代码如下：

%%图像的傅里叶变换%%
imA=imread('rice.tif','tif'); %读取图像
imB=imread('lena.tiff','tif');
subplot(2,3,1);
imshow(imA);
title('原图像A');
subplot(2,3,2);
imshow(imB);
title('原图像B');
FA=fft2(imA);%对图像进行傅里叶变换
FB=fft2(imB);
fA=fftshift(FA); %对图像频谱进行移动，是0频率点在中心
fB=fftshift(FB);
sA=log(abs(fA));%获得傅里叶变换的幅度谱
sB=log(abs(fB));
phA=log(angle(fA)*180/pi);%获得傅里叶变换的相位谱
phB=log(angle(fB)*180/pi);
subplot(2,3,3);
imshow(sA,[]); %显示图像的度谱，参数与[]是为了将sA的值线形拉伸
title('图像A的傅里叶变换幅度谱');
subplot(2,3,4);
imshow(phA,[]); %显示图像傅里叶变换的相位谱
title('图像A傅里叶变换的相位谱');
subplot(2,3,5);
imshow(sB,[])
title('图像B的傅里叶变换幅度谱');
subplot(2,3,6);
imshow(phB,[]);
title('图像B傅里叶变换的相位谱');
A=ifft2(FA);%傅里叶反变换
B=ifft2(FB);
figure
subplot(1,2,1);
imshow(A,[]);
title('傅里叶反变换得到的A图像');
subplot(1,2,2);
imshow(B,[]);
title('傅里叶反变换的到的B图像');

结果如下：