理论基础

时域：以时间为横坐标

频域：以频率的倒数为横坐标，可以看出，频域更加简单。



相位：与时间差有关的一个概念。
傅里叶说，任何连续周期信号，可以由一组适当的正弦曲线组合而成。我们知道，正弦曲线可以转换为频域信号，所以：任何连续周期信号，都可以转换成频域信号。并且这个过程是可逆的。



程序实现

1. 傅里叶变换
numpy.fft.fft2


实现傅里叶变换。
返回一个复数数组。





numpy.fft.fftshift
效果如图所示




将零频率分量移到频谱中心。



20*np.log( np.abs( fshift ) )

将傅里叶变换的计算结果映射到【0,255】这个区间内。
import cv2
import numpy as np
import matplotlib.pyplot as plt

o=cv2.imread('image\\equ2.bmp',0)
f=np.fft.fft2(o)                              #傅里叶变换
fshift=np.fft.fftshift(f)                     #零频率移到中心
result= 20 * np.log(np.abs(fshift))           #阈值转换
plt.subplot(121),plt.imshow(o,cmap='gray'),plt.title('original'),plt.axis('off')
plt.subplot(122),plt.imshow(result,cmap='gray'),plt.title('result'),plt.axis('off')
plt.show()




2. 逆傅里叶变换
numpy.fft.ifft2


实现逆傅里叶变换。
返回一个复数数组。

numpy.fft.ifftshift

fftshift的逆函数,将低频从中心移到左上角。
iimg=np.abs( 逆傅里叶变换结果)

设置值得范围
将图像进行傅里叶变换后，再进行逆傅里叶变换，与原图片对比。
import cv2
import numpy as np
import matplotlib.pyplot as plt

o=cv2.imread('image\\boat.bmp',0)
f=np.fft.fft2(o)
fshift=np.fft.fftshift(f)                     #傅里叶变换
ishift=np.fft.ifftshift(fshift)
io=np.fft.ifft2(ishift)
io=np.abs(io)                                 #逆傅里叶变换
plt.subplot(121),plt.imshow(o,cmap='gray'),plt.title('original'),plt.axis('off')
plt.subplot(122),plt.imshow(io,cmap='gray'),plt.title('result'),plt.axis('off')
plt.show()

/*
opencv的傅里叶变换及逆变换实现。傅里叶变换需要将数据表示为复数形式，通过一个两通道矩阵来记录复数的实部和虚部，然后通过cvDFT来实现变换。对于图片变换也是一样，只是矩阵换成IplImage 格式，用两个单通道图片来分别表示实部和虚部，用两通道图片来存放变换结果。
lbd.uestc@http://www.doczj.com/doc/f8c8d5c62cc58bd63186bda0.html
*/
#include
#include
#include
#include
usingnamespace std;
int main(int argc,char*argv[])
{
float r[] = {23,23,214,214,32,53};
float i[] = {22,32,33,42,31,12};
CvMat *re = cvCreateMat(1,6,CV_32FC1);//实部
CvMat *im = cvCreateMat(1,6,CV_32FC1);//虚部
CvMat *re_dst = cvCreateMat(1,6,CV_32FC1);//变换后实部
CvMat *im_dst = cvCreateMat(1,6,CV_32FC1);//变换后虚部
CvMat *re_ = cvCreateMat(1,6,CV_32FC1);//反变换后实部
CvMat *im_ = cvCreateMat(1,6,CV_32FC1);//反变换后虚部
cvInitMatHeader(re, 1, 6, CV_32FC1, r); //初始化实部
cvInitMatHeader(im, 1, 6, CV_32FC1, i); //初始化虚部
CvMat *sum_src = cvCreateMat(1,6,CV_32FC2);//实部虚部结合体
CvMat *sum_dst = cvCreateMat(1,6,CV_32FC2);//变换后实部虚部结合体
CvMat *sum_ = cvCreateMat(1,6,CV_32FC2);//反变换后实部虚部结合体
cout<
cout<
for(int j = 0; j < re->cols; j++)
{
cout << cvmGet(re,0,j)
cout<
for(int i = 0; i < im->cols; i++)

积分线性变换可以将信号或图像在更适合的域内表达，并且使得解决相关问题更容易，在图像分析中最常用的积分显示变换是傅里叶变换、离散余弦变换与小波变换。

1d傅里叶变换由傅里叶(Fourier)提出，1d傅里叶变换将函数f(x)变换到频率域F(t)表达。F称作频谱，可以显示不同频率的相对成分。用i(根号下-1)表示虚数单位，1d连续傅里叶变换表达为：



逆变换为



相应地，对于下标从0开始的离散变量f(x)，1d离散傅里叶变换和逆变换分别为：



离散傅里叶变换(DFT)的时间复杂度为O(n2)，使用快速傅里叶变换(FFT)会降低为 O(nlogn)

从1d傅里叶变换可容易推广到2d，对于连续空间f，2d连续傅里叶变换和逆变换如下：



(x,y)表示图像坐标，(i,j)表示空间频率。

离散傅里叶变换和逆变换如下：



在matlab中，可使用fft2进行2d快速变换(1d使用fft)，如果用基本公式实现，复杂度为O（n4）。为了验证效果一下裁剪了局部：

I=imread('lena.jpg');
I=I(1:50,1:50);
[m,n]=size(I);
F=zeros(m,n);
%% FT
for u=1:m
    for v=1:n
        for x=1:m
            for y=1:n
                F(u,v)=F(u,v)+double(I(x,y))*exp(-2*pi*i*(u*x/m+v*y/n));
            end
        end
    end
end
I2=zeros(m,n);
for x=1:m
    for y=1:n
        for u=1:m
            for v=1:n
                I2(x,y)=I2(x,y)+double(F(u,v))*exp(2*pi*i*(u*x/m+v*y/n))/(m*n);
            end
        end
    end
end
I2=uint8(I2);
subplot(121),imshow(I);
subplot(122),imshow(I2);
F2=fft2(I);%直接调用
figure,
F=abs(F);
F2=abs(F2);
subplot(121),imshow(F,[]);
subplot(122),imshow(F2,[]);
F=fftshift(F);
F2=fftshift(F2);
F=log(F+1);
F2=log(F2+1);
figure,
subplot(221),imshow(F,[]);
subplot(222),imshow(F2,[]);
subplot(223),mesh(F);
subplot(224),mesh(F2);


其中，I为原图像，F为I的频谱，F2为直接fft2的频率，I2为F的逆变换结果(也可通过ifft2得到)，i为系统默认设置的常数跟下-1。直接运算或fft得到的F（如figure2）的频率高点都分布在四个角，直接使用fftshift可以把四个角的高频信息移动到最中间(imshow(I,[])在显示图像时会自动调整归一化，让图像正确显示，或者手动归一化实现后inshow(I))。之后abs将负实数和虚数部分调整为正实数，由于高频低频差距过大，所以用log(F+1)缓和高低。运行结果分别如下：







 figure1为裁剪后的图像与进行变换和逆变换后的图像，figure2为通过代码实现和直接调用库得到的未调整前的频谱(左图和右图角度不同，但内容一样)，figure3为调整后的频谱以及他们的mesh面。

将图像由空间域转为频率域在图像压缩以及去噪等方面有很好的应用。

之后文章：DCT WT

 

代码实现了快速傅里叶变换及其逆变换的计算，幅值、初始相位的计算
import numpy as np
import matplotlib.pyplot as plt
Fs=500#采样频率
ts=1/Fs#采样时间间隔
L=500#信号长度
t=np.linspace(0,L-1,L)*ts

x=2+0.7*np.cos(2*np.pi*50*t+np.pi/2)+2*np.cos(2*np.pi*200*t+np.pi/6)
y=x#+2*np.random.rand()
fft_y=np.fft.fft(y)#做傅里叶变换
yf=np.abs(y)
plt.rcParams['font.sans-serif']=['SimHei']#中文乱码处理
plt.rcParams['axes.unicode_minus']=False
plt.subplot(511)
plt.plot(t[0:100],x[0:100])
plt.title("原始信号")
plt.legend(("采样频率Fs=1000"))
plt.subplot(512)
f=np.linspace(0,1,int(L))*Fs
k=np.arange(0,L)
ff=Fs*k/L-Fs/2
plt.plot(f[1:],2*np.abs(fft_y[1:])/L)
plt.grid(True)
xi=np.abs(np.fft.ifft(fft_y))
plt.subplot(513)
plt.plot(t[0:100],xi[0:100])
#plt.show()
#plt.title("")
#####下面计算直流分量,a1[0]就是直流分量
a1=abs(fft_y)/L#直流分量计算
plt.subplot(514)
plt.plot(a1)
#plt.show()
#####下面计算初始相位的值
a=np.real(fft_y)
b=np.imag(fft_y)
p=np.arctan2(b,a)
p=p*180/np.pi
plt.subplot(515)
plt.plot(f[1:int(L/2)],p[1:int(L/2)])
plt.show()