图像处理-傅里叶变换

傅里叶变换是一种函数在空间域和频率的变换，从空间域到频率域的变换是傅里叶变换，而从频率域到空间域的转换叫做傅里叶的反变换

时域和频域：

1、频域

是指对函数或信号进行分析时，分析其和频率有关的部分，而不是和时间有关的部分，和时域相对

2、时域

是描述数学函数或者物理信号对时间的关系。例如一个信号的时域波形可以表示信号随时间的变化，在研究时域的信号时，常用示波器将信号转换为其时域的波形

3、两者之间的关系

时域（信号对时间的函数）和频率（信号对频率的函数）的变换在数学上通过积分变化实现的，对周期信号可以直接使用傅里叶变换，对非周期信号则要进行周期扩展，使用拉普拉斯变换

信号在频率域中的表现：

在频域中，频率越大说明原始信号变化速度越快；频率越小表示原始信号越平缓；当频率为0时，表示直流信号，没有变换。故频率的大小反映了信号的变换快慢。高频分量解释信号的突变部分，而低频分量决定信号的整体形象。

在图像处理中，频率反映了图像在空域灰度变换剧烈程度，也就是图像灰度的变化速度，也就是图像的梯度大小。

对图像而言，图像的边缘部分是突变部分，变化较快，故反映在频域上是高频分量

图像的噪声大部分情况下是高频部分

图像平缓变化部分则为低频分量

即傅里叶变换提供一个角度来观察图像，可以将图像从灰度分布转换到频率分布上来观察图像的特征

重要概念：

1、图像高频分量

图像突变部分，在某些情况下指图像边缘信息，噪声，或者两者混合

2、低频分量

图像变换平缓的部分，即图像的轮廓

3、高通滤波器

让图像使低频分量抑制，高频通过

4、低通滤波器

与高通相反，让图像使高频分量抑制，低频分量通过

5、带通滤波器

使图像在某部分的频率通过，其他过低或过高都抑制

6、带阻滤波器

与带通相反

傅里叶变换

也称作傅立叶变换，表示能将满足一定条件的某个函数表示成三角函数（正弦和/或余弦函数）或者它们的积分的线性组合。在不同的研究领域，傅立叶变换具有多种不同的变体形式，如连续傅立叶变换和离散傅立叶变换。傅里叶变换是一种分析信号的方法，它可分析信号的成分，也可用这些成分合成信号。许多波形可作为信号的成分，比如正弦波、方波、锯齿波等，傅里叶变换用正弦波作为信号的成分。
傅里叶变换的实质是将一个信号分离为无穷多多正弦/复指数信号的加成，也就是说，把信号变成正弦信号相加的形式——既然是无穷多个信号相加，那对于非周期信号来说，每个信号的加权应该都是零——但有密度上的差别，你可以对比概率论中的概率密度来思考一下——落到每一个点的概率都是无限小，但这些无限小是有差别的所以，傅里叶变换之后，横坐标即为分离出的正弦信号的频率，纵坐标对应的是加权密度

作用

举例说明：傅里叶变换可以将一个时域信号转换成在不同频率下对应的振幅及相位，其频谱就是时域信号在频域下的表现，而反傅里叶变换可以将频谱再转换回时域的信号。

最简单最直接的应用就是时频域转换，比如在移动通信的LTE系统中，要把接收的信号从时域变成频域，就需要使用FFT（快速傅里叶变换）。

又例如对一个采集到的声音做傅立叶变化就能分出好几个频率的信号。比如南非世界杯时，南非人吹的呜呜主拉的声音太吵了，那么对现场的音频做傅立叶变化（当然是对声音的数据做），会得到一个展开式，然后找出呜呜主拉的特征频率，去掉展开式中的那个频率的sin函数，再还原数据，就得到了没有呜呜主拉的嗡嗡声的现场声音。

而对图片的数据做傅立叶，然后增大高频信号的系数就可以提高图像的对比度。

同样，相机自动对焦就是通过找图像的高频分量最大的时候，就是对好了

频域中图像增强

• 可以利用频率成分和图像外表之间的对应关系。一些在空间域表述困难的增强任务，在频率域中变得非常普通
• 滤波在频率域更为直观，它可以解释空间域滤波的某些性质
• 可以在频率域指定滤波器，做反变换，然后在空间域使用结果滤波器作为空间域滤波器的指导
• 一旦通过频率域试验选择了空间滤波，通常实施都在空间域进行

反傅里叶变换

一维连续傅里叶变换及反变换

 单变量连续函数f(x)的傅里叶变换F(u)定义为：

其中，j = $sqrt {{ m{ - }}1} $=±i
给定F(u),通过傅里叶反变换可以得到f(x)：

二维连续傅里叶变换及反变换

二维连续函数f(x,y)的傅里叶变换F(u,v)定义为：

如果f(x,y)是实函数，它的傅里叶变换是对称的，即　　F(u,v) = F(− u,−v)

傅里叶变换的频率谱是对称的　　|F(u,v)| =| F(− u,−v)|

给定F(u,v),通过傅里叶反变换可以得到 f(x,y)：

一维离散傅里叶变换及反变换

单变量离散函数f(x)(x=0,1,2,…,M-1)的傅里叶变换F(u)定义为：

其中，u=0,1,2,…,M-1
从欧拉公式：


给定F(u),通过傅里叶反变换可以得到f(x)：

其中，x=0,1,2,…,M-1

二维离散傅里叶变换及反变换

图像尺寸为M×N的函数f(x,y)的DFT（离散傅里叶变换）为：

u=0,1,2,…,M-1, v=0,1,2,…,N-1
给出F(u,v),可通过反DFT（离散傅里叶变换）得到f(x,y)：

x=0,1,2,…,M-1, y=0,1,2,…,N-1
注：u和v是频率变量，x和y是空间域图像变量
F(0,0)表示：

这说明：假设f(x,y)是一幅图像，在原点的傅里叶变换等于图像的平均灰度级（M*N是总的像素点，f(x,y)是(x,y)点的灰度值，将所有的像素点的灰度值求和然后除以总的个数即为平均灰度值）

傅里叶变换的一维极坐标表示

幅度或频率谱为：

R(u)和I(u)分别是F(u)的实部和虚部
相角或相位谱为：

功率谱为：

f(x)的离散表示：

F(u)的离散表示：

傅里叶变换的二维极坐标表示

二维DFT的极坐标表示：

幅度或频率谱为：

R(u,v)和I(u,v)分别是F(u,v)的实部和虚部
相角或相位谱为：

功率谱为：

F(u,v)的原点变换：

用(-1)^x+y ✖ f(x,y),将F(u,v)原点变换到率坐标下的(M/2，N/2)，它是M×N区域的中心
u=0,1,2,…,M-1, v=0,1,2,…,N-1

傅里叶变换的性质

平移性

以⇔表示函数和其傅里叶变换的对应性

注：u和v是频率变量，x和y是空间域图像变量
公式（1）表明将f(x,y)与一个指数项相乘就相当于把其变换后的频域中心f(u,v) 移动到新的位置 f(u-uo,v-v0)
公式（2）表明将F(u,v)与一个指数项相乘就相当于把其变换后的空域中心f(x,y) 移动到新的位置 f(x-x0,y-y0)
公式（2）表明对f(x,y)的平移不影响其傅里叶变换的幅值

当u₀=M/2且v₀=N/2,

带入（1）和（2），得到

分配律

傅里叶变换对加法满足分配律，但对乘法则不满足：

尺度变换（缩放）

给定2个标量a和b，可以证明对傅里叶变换下列2个公式成立：

旋转性

引入极坐标 x = r cosθ, y = rsinθ,u =ω cosϕ,v =ωsinϕ
将f(x,y)和F(u,v)转换为 f (r,θ ) 和F(ω,ϕ)。将它们带入傅里叶变换对得到：

f(x,y)旋转角度θ 0，F(u,v)也将转过相同的角度
F(u,v)旋转角度θ 0，f(x,y)也将转过相同的角度

周期性和共轭对称性

尽管F(u,v)对无穷多个u和v的值重复出现，但只需根据在任一个周期里的N个值就可以从F(u,v)得到f(x,y)
只需一个周期里的变换就可将F(u,v)在频域里完全确定
同样的结论对f(x,y)在空域也成立

如果f(x,y)是实函数，则它的傅里叶变换具有共轭对称性

其中，F*(u,v)为F(u,v)的复共轭（当两个复数实部相等,虚部互为相反数时,这两个复数叫做互为共轭复数）

分离性

F(x,v)是沿着f(x,y)的一行所进行的傅里叶变换。当x=0,1,…,M-1，沿着f(x,y)的所有行计算傅里叶变换。

二维傅里叶变换的全过程

 先通过沿输入图像的每一行计算一维变换
 再沿中间结果的每一列计算一维变换
 可以改变上述顺序，即先列后行
 上述相似的过程也可以计算二维傅里叶反变换

平均值

由二维傅里叶变换的定义：

所以，

上式说明：如果f(x,y)是一幅图像，在原点的傅里叶变换即等于图像的平均灰度级

卷积理论

卷积是空间域过滤和频率域过滤之间的纽带
大小为M×N的两个函数f(x,y)和h(x,y)的离散卷积：

卷积定理：

相关性理论

相关的重要应用在于匹配：确定是否有感兴趣的物体区域
 f(x,y)是原始图像
 h(x,y)作为感兴趣的物体或区域（模板）
 如果匹配，两个函数的相关值会在h找到f中相应点的位置上达到最大
大小为M×N的两个函数f(x,y)和h(x,y)的相关性定义为：

f* 表示f的复共轭。对于实函数，f*＝f
相关定理：

自相关理论：

FFT

采用快速傅里叶变换（FFT）算法能使计算机计算离散傅里叶变换所需要的乘法次数大为减少，特别是被变换的抽样点数N越多，FFT算法计算量的节省就越显著。
函数或信号可以透过一对数学的运算子在时域及频域之间转换。和傅里叶变换作用一样

为什么需要快速傅里叶变换？

人们想让计算机能处理信号 但由于信号都是连续的、无限的，计算机不能处理，于是就有了傅里叶级数、傅里叶变换，将信号由时域变到频域，把一个信号变为有很多个不同频率不同幅度的正弦信号组成，这样计算机就能处理了，但又由于傅里叶变换中要用到卷积计算，计算量很大，计算机也算不过来，于是就有了快速傅里叶变换，大大降低了运算量，使得让计算机处理信号成为可能。快速傅里叶变换是傅里叶变换的快速算法而已，主要是能减少运算量和存储开销，对于硬件实现特别有利。

• 对u的M个值中的每一个都需进行M次复数乘法(将f(x) ✖  e− j2πux / M )和M-1次加法，即复数乘法和加法的次数都正比于M2
• 快速傅里叶变换(FFT)则只需要Mlog^2M次运算
• FFT算法与原始变换算法的计算量之比是log^2M/M，如M=1024≈103,则原始变换算法需要106次计算，而FFT需 要104次计算，FFT与原始变换算法之比是1：100

只考虑一维的情况，根据傅里叶变换的分离性可知，二维傅里叶变换可由连续2次一维傅里叶变换得到

基本思想

FFT算法基于一个叫做逐次加倍的方法。通过推导将原始傅里叶转换成两个递推公式：

公式推导

假设M的形式是
M = 2ⁿ， n为正整数。因此，M可以表示为：M = 2K 。将M=2K带入上式：

特性：

• 一个M个点的变换，能够通过将原始表达式分成两个部分来计算
• 通过计算两个（M/2）个点的变换。得Feven(u)和 Fodd(u)
• 奇部与偶部之和得到F(u)的前(M/2)个值
• 奇部与偶部之差得到F(u)的后(M/2)个值。且不需要额外的变换计算

归纳快速傅立叶变换的思想

1）通过计算两个单点的DFT，来计算两个点的DFT

2）通过计算两个双点的DFT，来计算四个点的DFT，…，以此类推

3）对于任何N=2m的DFT的计算，通过计算两个N/2点的DFT，来计算N个点的DFT

举例

设：有函数f(x)，其N = 2³ = 8,有：{f(0),f(1),f(2),f(3),f(4),f(5),f(6),f(7)}
计算：
{F(0),F(1),F(2),F(3),F(4),F(5),F(6),F(7)}

解法：
首先分成奇偶两组，有：
{ f(0), f(2), f(4), f(6) }
{ f(1), f(3), f(5), f(7) }
为了利用递推特性，再分成两组，有：
{ f(0), f(4) }， { f(2), f(6) }
{ f(1), f(5) }， { f(3), f(7) }
对输入数据的排序可根据一个简单的位对换规则进行：
如用x表示f(x)的1个自变量值，那么它排序后对应的值可通过把x表示成二进制数并对换各位得到。例如N=2³,f(6)排序后为f(3),因为6＝1102而0112 ＝3
把输入数据进行了重新排序，则输出结果是正确的次序。反之不把输入数据进行排序，则输出结果需要重新排序才能得到正确的次序
地址的排序：——按位倒序规则
例如：N = 2³ = 8

2）计算顺序及地址增量：2ⁿ, n = 0,1,2…

傅里叶变换的物理意义

1、图像的频率是表征图像中灰度变化剧烈程度的指标，是灰度在平面空间上的梯度。（灰度变化得快频率就高，灰度变化得慢频率就低）

如：大面积的沙漠在图像中是一片灰度变化缓慢的区域，对应的频率值很低；而对于地表属性变换剧烈的边缘区域在图像中是一片灰度变化剧烈的区域，对应的频率值较高。

傅立叶变换在实际中有非常明显的物理意义，设f是一个能量有限的模拟信号，则其傅立叶变换就表示f的谱。从纯粹的数学意义上看，傅立叶变换是将一个函数转换为一系列周期函数来处理的。从物理效果看，傅立叶变换是将图像从空间域转换到频率域，其逆变换是将图像从频率域转换到空间域。换句话说，傅立叶变换的物理意义是将图像的灰度分布函数变换为图像的频率分布函数，傅立叶逆变换是将图像的频率分布函数变换为灰度分布函数

2、傅立叶变换以前，图像（未压缩的位图）是由对在连续空间（现实空间）上的采样得到一系列点的集合，我们习惯用一个二维矩阵表示空间上各点，则图像可由z=f(x,y)来表示。由于空间是三维的，图像是二维的，因此空间中物体在另一个维度上的关系就由梯度来表示，这样我们可以通过观察图像得知物体在三维空间中的对应关系。为什么要提梯度？因为实际上对图像进行二维傅立叶变换得到频谱图，就是图像梯度的分布图，当然频谱图上的各点与图像上各点并不存在一一对应的关系，即使在不移频的情况下也是没有。傅立叶频谱图上我们看到的明暗不一的亮点，实际上图像上某一点与邻域点灰度值差异的强弱，即梯度的大小，也即该点的频率的大小（差异/梯度越大，频率越高，能量越低，在频谱图上就越 暗。差异/梯度越小，频率越低，能量越高，在频谱图上就越 亮。换句话说，频率谱上越亮能量越高，频率越低，图像差异越小/平缓）。一般来讲，梯度大则该点的亮度强，否则该点亮度弱。频谱图，也叫功率图。

在经过频谱中心化（用(-1)^x+y乘以输入的图像函数

）后的频谱中，中间最亮的点是最低频率，属于直流分量（DC分量）(当频率为0时，表示直流信号，没有变化。在原点(u,v两个频率域变量均为零)的傅里叶变换即等于图像的平均灰度级,F(0，0)称做频 率谱的直流成分)。越往边外走，频率越高。所以，频谱图中的四个角和X,Y轴的尽头都是高频，如下图：

我们首先就可以看出，图像的能量分布，如果频谱图中暗的点数更多，那么实际图像是比较柔和的（因为各点与邻域差异都不大，梯度相对较小），反之，如果频谱图中亮的点数多，那么实际图像一定是尖锐的，边界分明且边界两边像素差异较大的。对频谱移频到原点以后，可以看出图像的频率分布是以原点为圆心，对称分布的。将频谱移频到圆心除了可以清晰地看出图像频率分布以外，还有一个好处，它可以分离出有周期性规律的干扰信号，比如正弦干扰，一副带有正弦干扰，移频到原点的频谱图上可以看出除了中心以外还存在以某一点为中心，对称分布的亮点集合，这个集合就是干扰噪音产生的，这时可以很直观的通过在该位置放置带阻滤波器消除干扰

实验

DFT数据加密

我们已知，利用复数来表达频谱

假设image size is 100*100

A = imread(image);
B = fft2(A); %此时，B的数据类型是复数，complex number
C = B + m*sth; %m是加权系数或者能量程度，sth 就是要加密的信息
D = ifft2(C);

上述四行代码做了个简单的加密工作，关键就是第三步，实际上是复数的相加，实际编码中，还可以考虑加在实部或者虚部，具体看要加密的信息格式和要求了。

同样，对应的解密工作

A = fft2(imread(image));
B = fft2(imread(D));
C = (B - A)/m; %复数相减
D = ifft2(C)

主要就是复数相减。

二维傅里叶fft2对简单矩阵的操作

二维快速傅里叶变换

    此 MATLAB 函数 使用快速傅里叶变换算法返回矩阵的二维傅里叶变换，这等同于计算 fft(fft(X).').'。如果 X 是一个多维数组，fft2
    将采用高于 2 的每个维度的二维变换。输出 Y 的大小与 X 相同。

    Y = fft2(X)
    Y = fft2(X,m,n)

%magic - 幻方矩阵

    此 MATLAB 函数 返回由 1 到 n2 的整数构成并且总行数和总列数相等的 n×n 矩阵。阶次 n 必须为大于或等于 3 的标量。

    M = magic(n)  ，n是维数
%

>> s=magic(2)

s =

     1     3
     4     2

>> f=fft2(s)

f =

    10     0
    -2    -4
%abs - 绝对值和复数的模

    此 MATLAB 函数 返回数组 X 中每个元素的绝对值。

    Y = abs(X)
%
>> a=abs(f)

a =

    10     0
     2     4

%3*3的矩阵的操作结果：

>> s=magic(3)

s =

     8     1     6
     3     5     7
     4     9     2

>> f=fft2(s)

f =

  45.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i
   0.0000 + 0.0000i  13.5000 + 7.7942i   0.0000 - 5.1962i
   0.0000 - 0.0000i   0.0000 + 5.1962i  13.5000 - 7.7942i

>> a=abs(f)

a =

   45.0000         0         0
    0.0000   15.5885    5.1962
    0.0000    5.1962   15.5885

%4*4的矩阵操作结果：

>> s=magic(4)

s =

    16     2     3    13
     5    11    10     8
     9     7     6    12
     4    14    15     1

>> f=fft2(s)

f =

   1.0e+02 *

   1.3600 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i
   0.0000 + 0.0000i   0.2000 + 0.0000i   0.0800 + 0.0800i   0.0000 - 0.1200i
   0.0000 + 0.0000i   0.3200 + 0.3200i   0.0000 + 0.0000i   0.3200 - 0.3200i
   0.0000 + 0.0000i   0.0000 + 0.1200i   0.0800 - 0.0800i   0.2000 + 0.0000i

>> a=abs(f)

a =

  136.0000         0         0         0
         0   20.0000   11.3137   12.0000
         0   45.2548         0   45.2548
         0   12.0000   11.3137   20.0000

二维快速傅里叶fft2对灰度图的操作

close all;
I=imread('input.jpg');%读入图像
J=rgb2gray(I);%将图像转换为灰度图
K=fft2(J);%对图像进行二维快速傅里叶变换
K=fftshift(K);%将频谱转移到中心，其实就是在傅里叶变换时乘以了某个因子
L=abs(K/256);%取模
%显示图片
figure;
subplot(131);
imshow(I);title('原图像')
subplot(132);
imshow(J);title('被转换为灰度图后')
subplot(133);
imshow(uint8(L));title('二维傅里叶变换后的频谱图')

I=imread('input.jpg');%读入图像
J=rgb2gray(I);%转换为灰度图
J=imrotate(J,45,'bilinear');	%将图像旋转45度角	
K=fft2(J);%对图像进行二维傅里叶比变换
K=fftshift(K);%转移频谱中心
L=abs(K/256);%取模
figure;
subplot(121);
imshow(J);title('原图像被旋转45度角')
subplot(122);
imshow(uint8(L));title('傅里叶变换后的图像')

I=imread('input.jpg');
J=rgb2gray(I);
J=imnoise(J, 'gaussian', 0, 0.01);%加入高斯噪声
K=fft2(J);
K=fftshift(K);
L=abs(K/256);
figure;
subplot(131);
imshow(I);title('原图像')
subplot(132);
imshow(J);title('高斯噪声污染')
subplot(133);
imshow(uint8(L));title('污染后进行傅里叶变换')

I=imread('input.jpg');
J=rgb2gray(I);
K=fft2(J);%傅里叶变换
L=fftshift(K);%转移频谱中心
L=abs(L/256);%取模
M=ifft2(K);%二维傅里叶反变换
figure;
subplot(131);
imshow(J);title('原灰度图像')
subplot(132);
imshow(uint8(L));title('傅里叶变换后的图像')
subplot(133);
imshow(uint8(M));title('反变换后的图像')

DFT嵌入水印

function p =watermark()
%使用DCT方法，实现图片水印嵌入
chars_info = char('123456789123456789123456789123456789');
fprintf('原始水印信息
 %s 
', chars_info);
bits_info = reshape(de2bi(uint8(chars_info),8,'left-msb')',[],1);
yuanImg=imread('Lena2.bmp');
k=50000;  %DFT用的系数

wmImg=[];
yuanImgYUV=[];
%彩色图像水印嵌入
yuanImgYUV = rgb2ycbcr(yuanImg);
%观察原始图像幅度频率特性
%fftAM=uint8(abs(fftshift(fft2(yuanImgYUV(:,:,1))))/100);
%imshow(fftAM)
yuanImgYUV(:,:,1) = EmbedWM_289(yuanImgYUV(:,:,1), bits_info, k,1);
wmImg = ycbcr2rgb(yuanImgYUV);
%观察嵌入水印后图像幅度频率特性
%fftAM=uint8(abs(fftshift(fft2(yuanImgYUV(:,:,1))))/100);
%figure(2);imshow(fftAM)
imwrite(wmImg,'wm.bmp');
figure;
subplot(1,2,1),imshow(yuanImg),title("原图");
subplot(1,2,2),imshow(wmImg),title("嵌入水印后");

function fw = EmbedWM_289(f, bits_info, k, type)
%本程序用来对图像进行DFT，通过改变DFT的分块系数的能量关系，实现水印的嵌入
%f:当前需要嵌入的图像
%bits_info:嵌入的水印内容
%k:嵌入水印的强度
%type：表示嵌入的类型，
%type=1时：(e1 - e) > k,表示水印0; (e2 - e1) > k,表示水印1
%type=2时：(e1 - e2) > k,表示水印1; (e2 - e1) > k,表示水印0
stepx = 4;%每bit嵌入区域的横向步长
stepy = 4;%每bit嵌入区域的纵向步长

 %k = 50000;

len_info = length(bits_info);%嵌入信息的长度

fftHuge = fftshift( fft2(f) );

if ( mod(size(fftHuge, 1), 2) == 0 )
    if ( mod(size(fftHuge, 2), 2) == 0 )
        fftHuge_ct = fftHuge(2:end, 2:end);
    else
        fftHuge_ct = fftHuge(2:end, :);
    end
else
    if ( mod(size(fftHuge, 2), 2) == 0 )
        fftHuge_ct = fftHuge(:, 2:end);
    else
        fftHuge_ct = fftHuge(:, :);
    end
end

dist_r = floor( size(fftHuge_ct, 1)/4 );
dist_c = floor( size(fftHuge_ct, 2)/4 );

fftHuge_wm = fftHuge_ct(dist_r+1: end-dist_r, dist_c+1: end-dist_c);

rmax_block = floor(size(fftHuge_wm,1)/(2*stepy));%纵向（行）最大分块数，由于中心对称性，故只处理上半部分
cmax_block = floor(size(fftHuge_wm,2)/stepx);%横向（列）最大分块数

if( rmax_block * cmax_block < len_info )
    disp('水印信息过长
');
    return;
end

SG1 = zeros(stepy*stepx/4, 2);
SG2 = zeros(stepy*stepx/4, 2);

% modified coefficients
% 右上角的系数
count = 1;
for r = 1 : stepy/2
    for c = stepx/2+1 : stepx
        SG1(count, 1) = r;
        SG1(count, 2) = c;
        count = count + 1;
    end
end

% 左下角的系数
count = 1;
for r = stepy/2+1 : stepy
    for c = 1 : stepx/2
        SG2(count, 1) = r;
        SG2(count, 2) = c;
        count = count + 1;
    end
end

count = 1;
for r = 1 : rmax_block
    for c = 1 : cmax_block
        fft_block = fftHuge_wm( (r-1)*stepy+1 : r*stepy, (c-1)*stepx+1 : c*stepx );
        %---------------------------系数修改-------------------------------- 
        % 计算能量及幅角
        ext = abs(fft_block); %幅度
        theta = angle(fft_block); %相角
        
        e1 = 0;
        e2 = 0;
        for i = 1 : size(SG1,1)
            e1 = e1 + ext(SG1(i,1), SG1(i,2));
            e2 = e2 + ext(SG2(i,1), SG2(i,2));
        end
        
        % 修改能量
        if(type==1)   %类型1
            if ( bits_info(count) == 0 && (e1 - e2) < k )
                delta = (k - e1 + e2)/(2*size(SG1,1));%每个系数的修改量，size(SG1,1)表示每个区域的点数
                for i = 1 : size(SG1,1)
                    ext(SG1(i,1), SG1(i,2)) = ext(SG1(i,1), SG1(i,2)) + delta;
                    ext(SG2(i,1), SG2(i,2)) = ext(SG2(i,1), SG2(i,2)) - delta;
                end 
            elseif ( bits_info(count) == 1 && (e2 - e1) < k)
                delta = (k - e2 + e1)/(2*size(SG1,1));%每个系数的修改量，size(SG1,1)表示每个区域的点数
                for i = 1 : size(SG1,1) 
                    ext(SG2(i,1), SG2(i,2)) = ext(SG2(i,1), SG2(i,2)) + delta;
                    ext(SG1(i,1), SG1(i,2)) = ext(SG1(i,1), SG1(i,2)) - delta;
                end
            end
        elseif(type==2)
            if ( bits_info(count) == 1 && (e1 - e2) < k )
                delta = (k - e1 + e2)/(2*size(SG1,1));%每个系数的修改量，size(SG1,1)表示每个区域的点数
                for i = 1 : size(SG1,1)
                    ext(SG1(i,1), SG1(i,2)) = ext(SG1(i,1), SG1(i,2)) + delta;
                    ext(SG2(i,1), SG2(i,2)) = ext(SG2(i,1), SG2(i,2)) - delta;
                end 
            elseif ( bits_info(count) == 0 && (e2 - e1) < k)
                delta = (k - e2 + e1)/(2*size(SG1,1));%每个系数的修改量，size(SG1,1)表示每个区域的点数
                for i = 1 : size(SG1,1) 
                    ext(SG2(i,1), SG2(i,2)) = ext(SG2(i,1), SG2(i,2)) + delta;
                    ext(SG1(i,1), SG1(i,2)) = ext(SG1(i,1), SG1(i,2)) - delta;
                end
            end
        end
        % 恢复FFT系数矩阵
        re = ext.*cos(theta);
        im = ext.*sin(theta);
        fft_block = re + 1i*im;
        %------------------------------------------------------------------
        fftHuge_wm( (r-1)*stepy+1 : r*stepy, (c-1)*stepx+1 : c*stepx ) = fft_block;%将修改后的FFT系数置回
        fftHuge_wm( end-r*stepy+1 : end-(r-1)*stepy, end-c*stepx+1 : end-(c-1)*stepx ) = rot90( conj(fft_block), 2 );%将修改后的共轭FFT系数转置后置回
        
        if ( count >= len_info )
            break;
        end
            
        count = count + 1;
    end
    
    if ( count >= len_info )
        break;
    end
end

fftHuge_ct(dist_r+1: end-dist_r, dist_c+1: end-dist_c) = fftHuge_wm;

if ( mod(size(fftHuge, 1), 2) == 0 )
    if ( mod(size(fftHuge, 2), 2) == 0 )
        fftHuge(2:end, 2:end) = fftHuge_ct;
    else
        fftHuge(2:end, :) = fftHuge_ct;
    end
else
    if ( mod(size(fftHuge, 2), 2) == 0 )
        fftHuge(:, 2:end) = fftHuge_ct;
    else
        fftHuge(:, :) = fftHuge_ct;
    end
end

fw = uint8( ifft2(ifftshift(fftHuge)) );

mse=mymse(f,fw,size(f,2),size(f,1));
psnr=10*log10(255^2/mse);
fprintf('PSNR = %f
', psnr);

输出：

原始水印信息
 123456789123456789123456789123456789 
PSNR = 43.801986

水印提取：

function p =extraWM()
%使用DCT方法，实现图片水印提取
%提取水印时不需要原始水印内容，但是需要水印的长度
chars_info = char('123456789123456789123456789123456789');
fprintf('原始水印信息
 %s 
', chars_info);
bits_info = reshape(de2bi(uint8(chars_info),8,'left-msb')',[],1);

%下面进行水印提取
wmImg=imread('wm.bmp');
%wmImg=imread('wm.jpg');
%wmImg=imread('wm2.bmp');
%彩色图像提取水印
wmImgYUV=[];
wmImgYUV = rgb2ycbcr(wmImg);
bits_info_ext = ExtractWM_289(wmImgYUV(:,:,1), length(bits_info),1);      
%二进制数组转换为字符串
extrastr='';
[m n]=size(bits_info_ext);
for x=1:m/8
    temp=bi2de(bits_info_ext((x-1)*8+1:x*8)','left-msb');
    extrastr=[extrastr char(temp)];
end
extrastr
ber(bits_info_ext, bits_info)

输出：

extrastr =

    '123456789123456789123456789123456789'


ans =

     0

参考

1、链接

2、数字图像处理(中译第三版)[冈萨雷斯]

3、链接

4、基于频域的信息加密----傅里叶变换