傅里葉變換是一種函數在空間域和頻率的變換，從空間域到頻率域的變換是傅里葉變換，而從頻率域到空間域的轉換叫做傅里葉的反變換
時域和頻域：
1、頻域
是指對函數或信號進行分析時，分析其和頻率有關的部分，而不是和時間有關的部分，和時域相對
2、時域
是描述數學函數或者物理信號對時間的關系。例如一個信號的時域波形可以表示信號隨時間的變化，在研究時域的信號時，常用示波器將信號轉換為其時域的波形
3、兩者之間的關系
時域（信號對時間的函數）和頻率（信號對頻率的函數）的變換在數學上通過積分變化實現的，對周期信號可以直接使用傅里葉變換，對非周期信號則要進行周期擴展，使用拉普拉斯變換
信號在頻率域中的表現：
在頻域中，頻率越大說明原始信號變化速度越快；頻率越小表示原始信號越平緩；當頻率為0時，表示直流信號，沒有變換。故頻率的大小反映了信號的變換快慢。高頻分量解釋信號的突變部分，而低頻分量決定信號的整體形象。
在圖像處理中，頻率反映了圖像在空域灰度變換劇烈程度，也就是圖像灰度的變化速度，也就是圖像的梯度大小。
對圖像而言，圖像的邊緣部分是突變部分，變化較快，故反映在頻域上是高頻分量
圖像的噪聲大部分情況下是高頻部分
圖像平緩變化部分則為低頻分量
即傅里葉變換提供一個角度來觀察圖像，可以將圖像從灰度分布轉換到頻率分布上來觀察圖像的特征
重要概念：
1、圖像高頻分量
圖像突變部分，在某些情況下指圖像邊緣信息，噪聲，或者兩者混合
2、低頻分量
圖像變換平緩的部分，即圖像的輪廓
3、高通濾波器
讓圖像使低頻分量抑制，高頻通過
4、低通濾波器
與高通相反，讓圖像使高頻分量抑制，低頻分量通過
5、帶通濾波器
使圖像在某部分的頻率通過，其他過低或過高都抑制
6、帶阻濾波器
與帶通相反
傅里葉變換
也稱作傅立葉變換，表示能將滿足一定條件的某個函數表示成三角函數（正弦和/或余弦函數）或者它們的積分的線性組合。在不同的研究領域，傅立葉變換具有多種不同的變體形式，如連續傅立葉變換和離散傅立葉變換。傅里葉變換是一種分析信號的方法，它可分析信號的成分，也可用這些成分合成信號。許多波形可作為信號的成分，比如正弦波、方波、鋸齒波等，傅里葉變換用正弦波作為信號的成分。
傅里葉變換的實質是將一個信號分離為無窮多多正弦/復指數信號的加成，也就是說，把信號變成正弦信號相加的形式——既然是無窮多個信號相加，那對於非周期信號來說，每個信號的加權應該都是零——但有密度上的差別，你可以對比概率論中的概率密度來思考一下——落到每一個點的概率都是無限小，但這些無限小是有差別的所以，傅里葉變換之后，橫坐標即為分離出的正弦信號的頻率，縱坐標對應的是加權密度
作用
舉例說明：傅里葉變換可以將一個時域信號轉換成在不同頻率下對應的振幅及相位，其頻譜就是時域信號在頻域下的表現，而反傅里葉變換可以將頻譜再轉換回時域的信號。
最簡單最直接的應用就是時頻域轉換，比如在移動通信的LTE系統中，要把接收的信號從時域變成頻域，就需要使用FFT（快速傅里葉變換）。
又例如對一個采集到的聲音做傅立葉變化就能分出好幾個頻率的信號。比如南非世界杯時，南非人吹的嗚嗚主拉的聲音太吵了，那么對現場的音頻做傅立葉變化（當然是對聲音的數據做），會得到一個展開式，然后找出嗚嗚主拉的特征頻率，去掉展開式中的那個頻率的sin函數，再還原數據，就得到了沒有嗚嗚主拉的嗡嗡聲的現場聲音。
而對圖片的數據做傅立葉，然后增大高頻信號的系數就可以提高圖像的對比度。
同樣，相機自動對焦就是通過找圖像的高頻分量最大的時候，就是對好了
頻域中圖像增強
- 可以利用頻率成分和圖像外表之間的對應關系。一些在空間域表述困難的增強任務，在頻率域中變得非常普通
- 濾波在頻率域更為直觀，它可以解釋空間域濾波的某些性質
- 可以在頻率域指定濾波器，做反變換，然后在空間域使用結果濾波器作為空間域濾波器的指導
- 一旦通過頻率域試驗選擇了空間濾波，通常實施都在空間域進行
反傅里葉變換
一維連續傅里葉變換及反變換
單變量連續函數f(x)的傅里葉變換F(u)定義為：
其中，j = $\sqrt {{\rm{ - }}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)的平移不影響其傅里葉變換的幅值
當u0=M/2且v0=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)則只需要Mlog2M次運算
- FFT算法與原始變換算法的計算量之比是log2M/M，如M=1024≈103,則原始變換算法需要106次計算，而FFT需 要104次計算，FFT與原始變換算法之比是1：100
只考慮一維的情況，根據傅里葉變換的分離性可知，二維傅里葉變換可由連續2次一維傅里葉變換得到
基本思想
FFT算法基於一個叫做逐次加倍的方法。通過推導將原始傅里葉轉換成兩個遞推公式：
公式推導
假設M的形式是
M = 2n， 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 = 23 = 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=23,f(6)排序后為f(3),因為6＝1102而0112 ＝3
把輸入數據進行了重新排序，則輸出結果是正確的次序。反之不把輸入數據進行排序，則輸出結果需要重新排序才能得到正確的次序
地址的排序：——按位倒序規則
例如：N = 23 = 8
2）計算順序及地址增量：2n, 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('原始水印信息\n %s \n', 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('水印信息過長\n');
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\n', psnr);
輸出：
原始水印信息 123456789123456789123456789123456789 PSNR = 43.801986
水印提取：
function p =extraWM()
%使用DCT方法，實現圖片水印提取
%提取水印時不需要原始水印內容，但是需要水印的長度
chars_info = char('123456789123456789123456789123456789');
fprintf('原始水印信息\n %s \n', 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、鏈接