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