1. 為什么頻域分析值得你花時間搞明白很多做圖像處理的朋友拿到一張圖第一反應(yīng)就是調(diào)亮度、拉對比度、套濾鏡這些都是在空間域里折騰像素。但如果你只會空間域那一套遇到周期性噪聲、摩爾紋、特定方向的干擾條紋基本就只能干瞪眼。傅里葉變換就是給你開了一扇窗——把圖像從像素排布的視角切換到頻率成分的視角很多在空間域里看起來無解的難題到了頻域就是一刀切的事。我第一次真正體會到傅里葉變換的威力是處理一批工業(yè)相機拍的金屬表面圖像。圖像上有一層規(guī)律性的細密條紋用中值濾波、高斯濾波怎么調(diào)都去不干凈要么殘留明顯要么把表面紋理也糊掉了。后來把圖轉(zhuǎn)到頻域一看那層條紋就是頻譜上幾個孤立的亮點用陷波濾波器精準(zhǔn)打掉逆變換回來干干凈凈紋理保留完好。從那以后我就養(yǎng)成了習(xí)慣空間域搞不定的先轉(zhuǎn)頻域看看。這篇文章面向的是已經(jīng)會用Matlab做基本圖像操作、但對頻域分析還比較陌生的讀者。我會從傅里葉變換在圖像處理中的實際意義講起把頻譜圖的讀法、濾波器的設(shè)計、代碼的逐行邏輯都拆開說清楚最后給出幾個我實際項目中反復(fù)用到的完整代碼模板。你不需要有深厚的數(shù)學(xué)功底但需要愿意動手把代碼跑一遍——頻域這東西光看是看不懂的必須自己調(diào)參數(shù)、看頻譜變化才能建立起直覺。2. 圖像傅里葉變換的物理意義與頻譜圖讀法2.1 從一維信號到二維圖像頻率到底代表什么一維信號的傅里葉變換大家應(yīng)該不陌生任何一段波形都可以分解成不同頻率正弦波的疊加。頻率高代表變化快頻率低代表變化慢。到了二維圖像這個概念需要擴展一下——圖像是二維信號所以頻率也是二維的分為水平方向的頻率和垂直方向的頻率。具體來說圖像中變化劇烈的區(qū)域?qū)?yīng)高頻邊緣、噪點、細紋理這些地方像素值在相鄰位置跳變很大。而變化平緩的區(qū)域?qū)?yīng)低頻大面積的天空、墻面、平滑漸變這些地方像素值在空間上緩慢變化。直流分量頻率為零就是整張圖的平均亮度。這里有個容易繞暈的點二維頻率是一個向量有方向和大小。水平條紋沿y方向變化劇烈對應(yīng)垂直方向的高頻垂直條紋對應(yīng)水平方向的高頻。斜條紋則對應(yīng)斜方向的高頻。理解這一點你才能看懂頻譜圖上那些亮點的位置含義。2.2 頻譜圖為什么要把零頻搬到中心Matlab的fft2函數(shù)輸出的頻譜零頻分量在矩陣的左上角。直接imshow出來會看到四個角是亮的中心是暗的這不符合人的直覺。所以標(biāo)準(zhǔn)操作是用fftshift把零頻搬到中心這樣頻譜圖中心代表低頻、越往外代表越高頻讀起來就直觀了。但這里有個坑fft2輸出的頻譜幅度范圍極大直流分量可能是其他分量的幾萬倍。直接顯示的話除了中心一個白點其他什么都看不見。所以必須做對數(shù)變換log(1 abs(F))。加1是為了避免log(0)的問題。這個操作叫對數(shù)幅度譜是頻域分析的標(biāo)準(zhǔn)可視化手段。我見過不少初學(xué)者在這一步卡住以為自己的頻譜算錯了其實就是沒做對數(shù)壓縮。記住頻譜可視化必須做對數(shù)變換否則你看到的永遠是一個白點加一片黑。2.3 頻譜圖的對稱性與實圖像的共軛對稱實值圖像的傅里葉變換具有共軛對稱性F(u,v) F*(-u,-v)。這意味著幅度譜關(guān)于中心點是中心對稱的。你在頻譜圖上看到的每一個亮點在中心對稱位置一定有一個同樣亮的點。這個性質(zhì)在濾波器設(shè)計時很有用——你只需要設(shè)計一半的濾波器另一半對稱過去就行能省不少事。另外這個對稱性也意味著頻譜圖的信息量實際上只有一半是獨立的。做壓縮的時候可以利用這一點不過那是另一個話題了。3. 頻域濾波的核心邏輯濾波器怎么設(shè)計才對3.1 理想低通、巴特沃斯、高斯三種低通濾波器的取舍低通濾波是頻域濾波里最基礎(chǔ)的操作目的是保留低頻、抑制高頻效果就是平滑去噪。但低通的具體實現(xiàn)方式不同效果差異很大。理想低通濾波器Ideal LPF的傳遞函數(shù)是一個硬截斷距離中心小于D0的頻率全部保留大于D0的全部置零。聽起來很完美但問題在于它的逆變換會產(chǎn)生振鈴效應(yīng) ringing。因為理想濾波器在頻域是矩形窗對應(yīng)空間域是sinc函數(shù)sinc函數(shù)的旁瓣就是振鈴的來源。圖像上表現(xiàn)為邊緣附近出現(xiàn)一圈圈波紋非常難看。巴特沃斯低通濾波器Butterworth LPF用一個平滑的過渡帶代替硬截斷階數(shù)n控制過渡帶的陡峭程度。n越大越接近理想濾波器振鈴越明顯n越小過渡越平滑振鈴越弱但高頻抑制不夠干凈。實際用的時候n取2到4之間比較平衡。高斯低通濾波器Gaussian LPF的傳遞函數(shù)是高斯函數(shù)它的逆變換也是高斯函數(shù)所以完全沒有振鈴。這是它最大的優(yōu)勢。代價是過渡帶比較寬高頻抑制的銳度不如巴特沃斯。如果你對振鈴特別敏感比如醫(yī)學(xué)圖像高斯是首選。下面這張表總結(jié)了三者的特性對比濾波器類型振鈴效應(yīng)過渡帶陡峭度計算復(fù)雜度適用場景理想低通嚴重極陡低理論分析實際少用巴特沃斯中等可調(diào)中通用去噪平衡選擇高斯無平緩低醫(yī)學(xué)圖像、對振鈴敏感場景3.2 高通濾波與邊緣提取的頻域視角高通濾波就是反過來抑制低頻、保留高頻。效果是突出邊緣和細節(jié)?？臻g域里的Sobel、Laplacian算子本質(zhì)上都是高通濾波只是實現(xiàn)方式不同。頻域做高通有個好處你可以精確控制保留哪個頻段的細節(jié)。比如你想提取中等尺度的紋理可以設(shè)計一個帶通濾波器只保留中間某段頻率。這是空間域卷積很難做到的。但高通濾波有個實際問題直流分量被抑制后圖像整體會變暗甚至出現(xiàn)負值。所以做完高通濾波后通常需要做歸一化或者加上一個偏移量。Matlab里可以用mat2gray自動歸一化到[0,1]范圍。3.3 陷波濾波精準(zhǔn)打擊周期性噪聲陷波濾波器Notch Filter是我個人用得最多的頻域工具。它的思路很簡單既然周期性噪聲在頻譜上表現(xiàn)為幾個孤立的亮點那我就針對這幾個點做抑制其他頻率不動。設(shè)計步驟是這樣的先看頻譜圖找到噪聲對應(yīng)的亮點位置記錄下它們的坐標(biāo)相對于中心點的偏移。然后構(gòu)造一個濾波器在這些位置附近設(shè)置一個抑制區(qū)域通常是圓形或矩形其他位置全為1。最后用這個濾波器乘以頻譜逆變換回去。這里的關(guān)鍵是確定抑制區(qū)域的半徑。半徑太小噪聲去不干凈半徑太大會誤傷周圍的正常頻率成分導(dǎo)致圖像出現(xiàn)新的偽影。我的經(jīng)驗是從較小的半徑開始試逐步增大直到噪聲剛好消失為止。通常半徑取5到15個像素之間。注意陷波濾波器的抑制區(qū)域必須關(guān)于中心點對稱。因為實圖像的頻譜是共軛對稱的你只抑制一個點逆變換回來會出現(xiàn)復(fù)數(shù)結(jié)果取實部后圖像會出問題。所以找到噪聲點后一定要同時抑制它的中心對稱點。4. 完整代碼實現(xiàn)與逐行拆解4.1 基礎(chǔ)傅里葉變換與頻譜可視化先把最基礎(chǔ)的流程跑通。下面這段代碼讀取一張圖像做傅里葉變換顯示頻譜圖。% 讀取圖像并轉(zhuǎn)為灰度 img imread(test_image.jpg); if size(img, 3) 3 img rgb2gray(img); end img im2double(img); % 二維傅里葉變換 F fft2(img); % 將零頻分量移到頻譜中心 F_shifted fftshift(F); % 計算幅度譜并做對數(shù)變換 magnitude abs(F_shifted); spectrum log(1 magnitude); % 顯示原圖和頻譜 figure; subplot(1, 2, 1); imshow(img); title(原始圖像); subplot(1, 2, 2); imshow(spectrum, []); title(對數(shù)幅度譜); colormap(gca, jet); colorbar;這段代碼里有幾個細節(jié)值得說。im2double把圖像轉(zhuǎn)成[0,1]的double類型這是必須的——如果直接用uint8做fft2Matlab會報錯或者結(jié)果不對。fftshift把零頻搬到中心前面已經(jīng)解釋過了。imshow(spectrum, [])里的空方括號是讓Matlab自動縮放顯示范圍因為對數(shù)譜的數(shù)值范圍不一定是[0,1]。跑完這段代碼你會看到頻譜圖中心有一個很亮的點直流分量然后向四周輻射出一些亮線。這些亮線的方向?qū)?yīng)圖像中邊緣的方向。如果圖像里有明顯的水平邊緣頻譜上就會有垂直方向的亮線反之亦然。4.2 高斯低通濾波的完整實現(xiàn)下面這段代碼實現(xiàn)高斯低通濾波包含濾波器構(gòu)造、濾波、逆變換的完整流程。% 假設(shè)img已經(jīng)讀取并轉(zhuǎn)為double類型 [M, N] size(img); % 傅里葉變換并中心化 F fft2(img); F_shifted fftshift(F); % 構(gòu)造高斯低通濾波器 D0 30; % 截止頻率需要根據(jù)圖像尺寸調(diào)整 u 0:M-1; v 0:N-1; [V, U] meshgrid(v, u); D sqrt((U - M/2).^2 (V - N/2).^2); H exp(-(D.^2) / (2 * D0^2)); % 濾波 G_shifted F_shifted .* H; % 逆變換 G ifftshift(G_shifted); img_filtered real(ifft2(G)); % 顯示結(jié)果 figure; subplot(1, 3, 1); imshow(img); title(原始圖像); subplot(1, 3, 2); imshow(H, []); title(高斯低通濾波器); subplot(1, 3, 3); imshow(img_filtered, []); title(濾波后圖像);這里的關(guān)鍵參數(shù)是D0也就是截止頻率。D0越小保留的低頻越少圖像越模糊D0越大保留的細節(jié)越多。對于512x512的圖像D0取30到60之間通常能看出明顯的平滑效果但不至于糊成一團。對于1024x1024的圖像D0要相應(yīng)放大到60到120。構(gòu)造濾波器時用meshgrid生成坐標(biāo)網(wǎng)格然后計算每個點到中心點的距離D。exp(-(D.^2)/(2*D0^2))就是高斯函數(shù)的標(biāo)準(zhǔn)形式。注意這里用的是D.^2而不是D因為高斯濾波器在頻域的傳遞函數(shù)是距離平方的函數(shù)。逆變換時ifftshift是fftshift的逆操作把零頻搬回左上角然后ifft2做逆變換。最后取real是因為浮點運算可能產(chǎn)生極小的虛部雖然理論上應(yīng)該是實數(shù)。4.3 陷波濾波去除周期性噪聲的實戰(zhàn)代碼這是我最常用的模板處理周期性噪聲非常有效。% 讀取含噪聲圖像 img imread(noisy_image.jpg); img im2double(img); [M, N] size(img); % 傅里葉變換 F fft2(img); F_shifted fftshift(F); % 顯示頻譜手動找噪聲點 figure; imshow(log(1 abs(F_shifted)), []); title(請觀察頻譜中的異常亮點); colormap(gca, jet); colorbar; % 假設(shè)通過觀察噪聲點位于以下位置相對于中心 % 這里需要你根據(jù)實際頻譜圖手動填寫 noise_points [45, 30; -45, -30; 30, -45; -30, 45]; % 示例坐標(biāo) % 構(gòu)造陷波濾波器 H ones(M, N); radius 10; % 抑制半徑 [V, U] meshgrid(1:N, 1:M); center_x N/2 0.5; center_y M/2 0.5; for k 1:size(noise_points, 1) % 噪聲點相對于中心的偏移 dx noise_points(k, 1); dy noise_points(k, 2); % 計算每個像素到噪聲點的距離 dist sqrt((U - (center_x dx)).^2 (V - (center_y dy)).^2); % 在噪聲點附近設(shè)置抑制 H(dist radius) 0; end % 濾波 G_shifted F_shifted .* H; G ifftshift(G_shifted); img_filtered real(ifft2(G)); % 顯示結(jié)果 figure; subplot(1, 2, 1); imshow(img); title(含噪聲圖像); subplot(1, 2, 2); imshow(img_filtered, []); title(陷波濾波后);這段代碼的核心是noise_points數(shù)組里面存儲了噪聲點相對于頻譜中心的坐標(biāo)偏移。你需要先跑一次頻譜顯示用鼠標(biāo)或者ginput函數(shù)找到亮點的位置然后填進去。radius控制抑制區(qū)域的大小從5開始試不夠再加。有個實用技巧如果噪聲點很多或者分布復(fù)雜可以寫一個交互式的腳本用ginput讓用戶點擊頻譜圖上的噪聲點自動記錄坐標(biāo)并生成濾波器。這樣效率高很多。5. 實操中容易踩的坑與排查思路5.1 圖像尺寸與fft2的性能陷阱fft2的計算復(fù)雜度跟圖像尺寸有關(guān)。如果圖像尺寸是2的冪次如256、512、1024計算速度最快。如果尺寸是質(zhì)數(shù)或者有大質(zhì)因子速度會慢很多。Matlab的fft2內(nèi)部會自動做padding優(yōu)化但如果你手動指定了變換長度最好選2的冪次。另一個坑是對彩色圖像做傅里葉變換時不能直接對三維數(shù)組做fft2。fft2只處理二維矩陣。正確做法是分別對R、G、B三個通道做變換處理完再合并?；蛘呦绒D(zhuǎn)成灰度圖如果顏色信息不重要的話。5.2 頻譜顯示全黑或全白的排查這是初學(xué)者最常遇到的問題。頻譜顯示全黑通常是因為沒做對數(shù)變換直流分量太大把其他分量壓沒了。頻譜顯示全白可能是數(shù)據(jù)類型問題——如果圖像是uint8abs之后數(shù)值可能超出顯示范圍。排查步驟先檢查class(img)是不是double再檢查有沒有做log(1abs(...))最后檢查imshow的第二個參數(shù)是不是空方括號。這三步走完99%的顯示問題都能解決。5.3 逆變換后圖像出現(xiàn)虛部或數(shù)值異常ifft2的結(jié)果理論上應(yīng)該是實數(shù)但浮點運算會引入極小的虛部。取real就行。但如果虛部很大說明濾波器設(shè)計有問題——最常見的原因是濾波器不滿足共軛對稱性。檢查你的濾波器是否關(guān)于中心點對稱特別是陷波濾波器有沒有漏掉對稱點。另一個可能的原因是濾波器值超出了[0,1]范圍。比如你用了H 1 - H_lowpass做高通結(jié)果可能有負值。負值濾波器在頻域是合法的但逆變換后圖像會出現(xiàn)負像素顯示時會被截斷。這時候需要做歸一化img_filtered mat2gray(img_filtered)。5.4 振鈴效應(yīng)的識別與抑制振鈴效應(yīng)表現(xiàn)為圖像邊緣附近出現(xiàn)波紋狀偽影。如果你用了理想低通濾波器振鈴幾乎必然出現(xiàn)。識別方法很簡單看邊緣附近有沒有周期性的明暗條紋。抑制振鈴的方法有三種一是換用高斯濾波器從根本上避免振鈴二是降低巴特沃斯濾波器的階數(shù)讓過渡帶更平滑三是在空間域做后處理比如用中值濾波輕微平滑一下。我通常首選高斯濾波器省心。6. 頻域分析還能怎么用幾個進階方向6.1 頻域做圖像配準(zhǔn)與相位相關(guān)傅里葉變換的相位信息包含了圖像的結(jié)構(gòu)信息。兩張有平移關(guān)系的圖像它們的傅里葉變換幅度譜相同相位譜相差一個線性項。通過計算互功率譜的逆變換可以找到峰值位置從而確定平移量。這就是相位相關(guān)法做圖像配準(zhǔn)的原理比空間域的相關(guān)法快很多而且對亮度變化不敏感。Matlab實現(xiàn)大概是這樣R ifft2(fft2(img1) .* conj(fft2(img2)))然后找R的峰值位置。峰值坐標(biāo)就是平移量。這個方法在遙感圖像配準(zhǔn)、醫(yī)學(xué)圖像對齊里用得很多。6.2 頻域濾波與卷積定理的等價性卷積定理說空間域的卷積等于頻域的乘法。所以你在空間域用卷積核做濾波等價于在頻域乘以對應(yīng)的傳遞函數(shù)。反過來也一樣。這意味著你可以自由選擇在哪個域操作——空間域卷積核小的時候直接卷積更快卷積核大的時候轉(zhuǎn)到頻域做乘法更高效。理解這個等價性對調(diào)試很有幫助。如果你在頻域設(shè)計的濾波器效果不對可以把它逆變換回空間域看看對應(yīng)的卷積核長什么樣往往能發(fā)現(xiàn)問題。6.3 小波變換與傅里葉變換的互補關(guān)系傅里葉變換給出的是全局頻率信息但丟失了空間位置信息——你只知道圖里有某個頻率成分但不知道它在哪個位置。小波變換通過多尺度分析同時保留了頻率和位置信息。對于非平穩(wěn)信號頻率成分隨時間或空間變化的信號小波變換往往比傅里葉變換更合適。不過小波變換的復(fù)雜度比傅里葉高不少Matlab的Wavelet Toolbox提供了完整的函數(shù)。如果你的應(yīng)用需要局部頻率分析可以了解一下。但大多數(shù)圖像濾波任務(wù)傅里葉變換已經(jīng)夠用了。7. 我在實際項目中的幾點體會頻域濾波這個工具入門容易精通難。我踩過的坑里最常見的就是濾波器參數(shù)調(diào)不好。D0、階數(shù)、抑制半徑這些參數(shù)沒有萬能值必須根據(jù)具體圖像反復(fù)試。我的習(xí)慣是寫一個帶滑動條的交互腳本邊調(diào)邊看效果找到最佳參數(shù)后再固化到代碼里。另一個體會是頻域濾波不是萬能的。它能精準(zhǔn)處理周期性噪聲和特定頻率的干擾但對于隨機噪聲比如高斯白噪聲頻域濾波的效果跟空間域平滑差不多沒有明顯優(yōu)勢。選工具要看問題類型別為了用而用。最后說一個實用技巧做頻域濾波之前先對圖像做一次頻譜分析看看能量分布。如果能量集中在低頻說明圖像比較平滑如果高頻能量也很多說明細節(jié)豐富或者噪聲大。這個分析能幫你判斷濾波器的截止頻率該設(shè)在哪里比盲目試參數(shù)高效得多。代碼方面我建議你把常用的濾波器構(gòu)造封裝成函數(shù)比如gaussianLPF(M, N, D0)、butterworthLPF(M, N, D0, n)、notchFilter(M, N, points, radius)。這樣每次用的時候直接調(diào)用不用重復(fù)寫坐標(biāo)計算。我自己的工具箱里就攢了這么一套函數(shù)用了好幾年省了很多時間。