1. 項(xiàng)目整體思路拆解一次搞定FFT去噪和數(shù)值積分先把這個(gè)標(biāo)題拆開(kāi)看里面其實(shí)是兩個(gè)完全獨(dú)立的經(jīng)典小項(xiàng)目一個(gè)是圖像去噪利用傅里葉變換把噪聲在頻域里“摳掉”另一個(gè)是求圓周率利用定積分把半圓面積算出來(lái)。它們唯一的共同點(diǎn)是都依賴Scipy這個(gè)庫(kù)一個(gè)用fft模塊一個(gè)用integrate模塊。這種組合式Demo特別適合正在學(xué)Python科學(xué)計(jì)算的人因?yàn)閮蓚€(gè)任務(wù)都足夠短小但背后牽扯的原理一個(gè)比一個(gè)深。先說(shuō)圖像去噪。傳統(tǒng)空間域處理比如中值濾波、均值濾波都是直接滑動(dòng)窗口去改像素值。傅里葉變換的思路完全不一樣它先把整張圖從“像素空間”變到“頻率空間”噪聲在頻率空間里通常表現(xiàn)為密密麻麻的高頻分量而圖像主體往往是低頻分量。我們?cè)陬l域里把高頻部分壓掉再逆變換回空間域噪聲就沒(méi)了。這個(gè)思路一旦打開(kāi)你會(huì)發(fā)現(xiàn)它不只是去噪還能用來(lái)做邊緣提取、圖像壓縮甚至水印去除。再說(shuō)定積分求圓周率。圓周率π是圓的周長(zhǎng)和直徑之比看起來(lái)和積分八竿子打不著。但如果你把單位圓的上半部分看成函數(shù)y sqrt(1 - x2)從-1到1對(duì)這個(gè)函數(shù)做定積分算出來(lái)的正好是半個(gè)圓的面積也就是π/2。那么π就等于這個(gè)積分結(jié)果乘以2。這個(gè)方案最妙的地方在于它把幾何問(wèn)題轉(zhuǎn)化成了數(shù)值計(jì)算問(wèn)題讓我們可以順理成章地引入scipy.integrate.quad這個(gè)自適應(yīng)積分器。兩個(gè)任務(wù)放在一起還有個(gè)額外好處FFT本身是一種精確變換誤差可控而數(shù)值積分天生就是近似計(jì)算誤差可以量化。這一精一近似形成對(duì)比能讓初學(xué)者更直觀地理解“數(shù)值計(jì)算到底在干什么”。我建議你按順序先做去噪再做積分因?yàn)槿ピ氲拇a更短視覺(jué)反饋也更直觀做完會(huì)特別有成就感。1.1 核心需求解析標(biāo)題背后到底想考你什么這個(gè)標(biāo)題看著像兩個(gè)小練習(xí)拼在一起其實(shí)考察點(diǎn)非常明確。第一部分考的是你是否理解傅里葉變換的物理意義而不是只會(huì)調(diào)用fft2然后瞎調(diào)參數(shù)。第二部分考的是你對(duì)定積分基本概念和自適應(yīng)性數(shù)值算法的掌握程度。如果你只是把代碼跑通那你最多只完成了30%剩下70%在于你能不能解釋清楚為什么噪聲在高頻為什么quad能算這種帶端點(diǎn)的開(kāi)根號(hào)函數(shù)為什么返回值里還帶著一個(gè)誤差估計(jì)從實(shí)操角度看還有一個(gè)隱含考察點(diǎn)環(huán)境搭建。Scipy不是Python自帶庫(kù)需要安裝numpy、scipy、matplotlib三個(gè)基礎(chǔ)包。很多人卡在去噪顯示圖像這一步其實(shí)只是沒(méi)裝matplotlib。我建議直接用pip install numpy scipy matplotlib一把梭然后用Jupyter Notebook來(lái)做因?yàn)槊恳徊蕉伎梢约磿r(shí)看到頻譜圖和濾波結(jié)果比寫(xiě)一個(gè)完整腳本調(diào)試要舒服太多。1.2 為什么是Scipy選型背后的真實(shí)原因有人可能會(huì)問(wèn)FFT用numpy.fft不就行了嗎為什么非要Scipy答案要從兩方面看。numpy.fft確實(shí)提供了基礎(chǔ)的fft2和ifft2夠用但接口比較糙而scipy.fft設(shè)計(jì)得更現(xiàn)代它對(duì)FFT的封裝支持worker線程、計(jì)劃緩存并且和老代碼的兼容性更好尤其在處理大規(guī)模數(shù)組時(shí)速度優(yōu)勢(shì)明顯。scipy.fft還有一個(gè)加分項(xiàng)它包含rfft2、irfft2這種專門(mén)處理實(shí)數(shù)輸入的快速路徑圖像恰恰是實(shí)數(shù)的用實(shí)數(shù)FFT能省下一半內(nèi)存和可觀的計(jì)算時(shí)間。積分這邊就更不用說(shuō)了。scipy.integrate.quad包裝的是經(jīng)典的數(shù)值積分庫(kù)內(nèi)部采用自適應(yīng)Gauss-Kronrod方法。你不用關(guān)心它內(nèi)部怎么選點(diǎn)它自己會(huì)在函數(shù)變化劇烈的地方加密采樣在平緩的地方減少采樣最終返回一個(gè)高精度的近似值和誤差估計(jì)。你要是自己寫(xiě)個(gè)矩形法或者辛普森法想要達(dá)到quad的精度計(jì)算量可能要差出好幾個(gè)數(shù)量級(jí)。我給一個(gè)非常直觀的建議先跟著本文用scipy.fft和scipy.integrate.quad跑通流程之后再去試試numpy.fft和手寫(xiě)積分對(duì)比一下性能和代碼復(fù)雜度你會(huì)對(duì)“為什么科學(xué)計(jì)算界普遍使用成熟庫(kù)”有切身體會(huì)。2. 傅里葉變換去噪從數(shù)學(xué)到代碼的完整落地2.1 圖像傅里葉變換的核心直覺(jué)想搞懂圖像FFT去噪先建立直覺(jué)別急著看公式。把圖像想象成一張起伏不平的地形圖像素亮度就是海拔。傅里葉變換做的事情就是把這整片地形拆解成無(wú)數(shù)不同頻率、不同方向的正弦波疊加。低頻正弦波是那種變化緩慢的大尺度起伏對(duì)應(yīng)圖像里漸變區(qū)域和主體輪廓高頻正弦波是那種變化劇烈的細(xì)小毛刺對(duì)應(yīng)圖像邊緣、紋理還有噪聲。噪聲為什么基本都在高頻因?yàn)樵肼暤奶攸c(diǎn)是隨機(jī)、突兀、局部突變比如椒鹽噪聲里那些黑白點(diǎn)和周?chē)袼亓炼炔町悩O大。從波的角度看這種突變只能靠高頻率的正弦波才能擬合出來(lái)。所以去噪的基本策略也就清楚了在頻域里把高頻成分衰減掉保留低頻主體再逆變換回去。這里有一個(gè)很多人忽略的點(diǎn)圖像邊緣本身也是高頻。低通濾波去噪的同時(shí)必然會(huì)讓邊緣變模糊。這是個(gè)物理上無(wú)法避免的權(quán)衡你能做的只是控制衰減強(qiáng)度在去噪和保留細(xì)節(jié)之間找一個(gè)平衡點(diǎn)。實(shí)操中我建議你每次改變?yōu)V波半徑后都輸出一張濾波后圖像觀察邊緣變化而不是只看一個(gè)PSNR數(shù)值。2.2 生成一張帶噪聲的測(cè)試圖先造一張干凈的測(cè)試圖再加噪聲這樣才能對(duì)比去噪前后的效果。我一般用圓盤(pán)圖原因是它有明顯的邊緣和平坦區(qū)域去噪效果肉眼一看便知。下面這段代碼生成一張256×256的圓盤(pán)圖然后加上椒鹽噪聲。import numpy as np import matplotlib.pyplot as plt from scipy import fft # 生成一個(gè)圓盤(pán)圖像 x np.linspace(-2, 2, 256) X, Y np.meshgrid(x, x) img 255 * ((X ** 2 Y ** 2) 1).astype(float) # 添加椒鹽噪聲 rng np.random.default_rng(42) noise rng.random(img.shape) img_noisy img.copy() img_noisy[noise 0.01] 0 # 隨機(jī)黑點(diǎn) img_noisy[noise 0.99] 255 # 隨機(jī)白點(diǎn) plt.figure(figsize(10, 5)) plt.subplot(1, 2, 1) plt.imshow(img, cmapgray) plt.title(Original) plt.axis(off) plt.subplot(1, 2, 2) plt.imshow(img_noisy, cmapgray) plt.title(Noisy) plt.axis(off) plt.show()我沒(méi)有用高斯噪聲而選椒鹽噪聲是因?yàn)榻符}噪聲在頻域里的特征更典型稀疏但強(qiáng)烈的沖激頻譜上呈現(xiàn)為散布在整個(gè)高頻區(qū)域的亮點(diǎn)。高斯噪聲則是均勻地鋪在所有頻率上處理起來(lái)更麻煩。對(duì)教學(xué)來(lái)說(shuō)椒鹽噪聲的濾波效果更明顯能讓你一眼看到魔法發(fā)生。2.3 頻域?yàn)V波核心代碼每個(gè)參數(shù)都有講究傅里葉變換去噪的標(biāo)準(zhǔn)流程五步二維FFT、頻譜中心化、構(gòu)造掩膜、頻域?yàn)V波、逆變換。每一步都有講究我拆開(kāi)說(shuō)。# 1. 二維FFT spectrum fft.fft2(img_noisy) # 2. 頻譜中心化把低頻移到圖像中心 spectrum_shifted fft.fftshift(spectrum) # 3. 構(gòu)造圓形低通掩膜 rows, cols img.shape crow, ccol rows // 2, cols // 2 radius 40 yy, xx np.ogrid[:rows, :cols] mask ((xx - ccol) ** 2 (yy - crow) ** 2) radius ** 2 # 4. 頻域?yàn)V波 filtered_spectrum spectrum_shifted * mask # 5. 逆變換取實(shí)部 filtered_img fft.ifft2(fft.ifftshift(filtered_spectrum)).real plt.figure(figsize(12, 4)) plt.subplot(1, 3, 1) plt.imshow(np.log(1 np.abs(spectrum_shifted)), cmapgray) plt.title(Spectrum) plt.axis(off) plt.subplot(1, 3, 2) plt.imshow(mask, cmapgray) plt.title(Mask (radius40)) plt.axis(off) plt.subplot(1, 3, 3) plt.imshow(filtered_img, cmapgray) plt.title(Filtered) plt.axis(off) plt.show()為什么要fftshift因?yàn)閒ft2算出來(lái)的頻譜低頻在四個(gè)角落高頻在中心附近這不符合我們“中心放低頻”的視覺(jué)習(xí)慣。fftshift把頻譜四角的低頻挪到中心方便畫(huà)圓掩膜。做完掩膜乘法后必須用ifftshift把頻譜還原回原來(lái)的排列方式再做逆變換這一步順序反了會(huì)得到一張莫名其妙的花屏。為什么要圓形掩膜而不是方形因?yàn)榉叫窝谀?huì)在濾波后產(chǎn)生明顯的振鈴效應(yīng)圖像邊緣出現(xiàn)一圈圈水波紋。圓形掩膜對(duì)應(yīng)的頻域邊界是連續(xù)旋轉(zhuǎn)對(duì)稱的振鈴會(huì)小得多。這是圖像處理里的老經(jīng)驗(yàn)新手經(jīng)常在這里翻車(chē)。radius40這個(gè)值不是拍腦袋定的它對(duì)應(yīng)的是保留頻率范圍。圖像尺寸256中心化后最大頻率半徑約為128。圓盤(pán)圖的主體是低頻半徑40已經(jīng)能覆蓋絕大多數(shù)主體頻率了。你可以試試調(diào)成20看效果圓盤(pán)邊緣會(huì)變糊但噪聲幾乎全沒(méi)了再調(diào)成80邊緣清晰了但噪聲點(diǎn)會(huì)殘留不少。這就是上一節(jié)說(shuō)的權(quán)衡。2.4 去噪效果評(píng)估別只顧著看漂不漂亮最關(guān)鍵的一步是評(píng)估效果。我常用的方法是同時(shí)輸出三張圖原圖、去噪圖、噪聲殘差圖。殘差就是原圖減去去噪圖理想情況下殘差里只能看到噪聲和一些邊緣細(xì)節(jié)如果殘差里出現(xiàn)明顯的圓形輪廓說(shuō)明低通濾波把邊緣信息也抹掉了。定量評(píng)估可以用信噪比SNR公式是信號(hào)能量除以噪聲能量。不過(guò)對(duì)于教學(xué)項(xiàng)目肉眼判斷反而更直接。我最喜歡的一個(gè)做法是放大局部區(qū)域?qū)Ρ缺热缬胮lt.xlim/ylim截取圓盤(pán)邊界附近放大能很清楚看到噪聲點(diǎn)是否被抹掉、邊緣是否變虛。還有一個(gè)小技巧低通濾波可以再做一次高通實(shí)驗(yàn)來(lái)加深理解。把掩膜取反只保留高頻再逆變換你會(huì)發(fā)現(xiàn)得到的圖幾乎是純邊緣和噪聲的疊加。這就能直觀印證“細(xì)節(jié)和噪聲都在高頻”這個(gè)結(jié)論。我強(qiáng)烈建議你跑一下這個(gè)對(duì)比實(shí)驗(yàn)比你單獨(dú)做低通濾波理解要深得多。3. 定積分求圓周率quad函數(shù)背后的數(shù)值藝術(shù)3.1 數(shù)學(xué)原理半圓面積就是π/2圓周率與定積分的關(guān)系本質(zhì)上是圓的面積公式。單位圓的方程是x2 y2 1取上半部分得到y(tǒng) sqrt(1 - x2)。對(duì)這個(gè)函數(shù)在[-1, 1]上做定積分得到的是上半圓的面積也就是π/2。所以π 2 × ∫_{-1}^{1} sqrt(1 - x2) dx。如果你覺(jué)得從-1到1麻煩也可以只算第一象限的部分x從0到1積分得到四分之一圓面積π/4最后乘4。兩種寫(xiě)法數(shù)值上沒(méi)差別但四分之一圓的寫(xiě)法積分區(qū)間更短對(duì)數(shù)值算法來(lái)說(shuō)通常更友好。這個(gè)題目看起來(lái)有點(diǎn)“繞”但它的教學(xué)價(jià)值極高它讓你意識(shí)到很多看起來(lái)和微積分無(wú)關(guān)的常數(shù)都可以通過(guò)數(shù)值積分算出來(lái)。你還可以擴(kuò)展出很多變體比如用同樣的方法算橢圓面積、球體體積只要寫(xiě)出被積函數(shù)quad就能處理。3.2 quad函數(shù)實(shí)操一行代碼算出半圓面積直接上代碼這是整個(gè)項(xiàng)目里最短的一段但信息量一點(diǎn)不少。from scipy.integrate import quad import numpy as np def half_circle(x): return np.sqrt(1 - x ** 2) area_half, error_estimate quad(half_circle, -1, 1) pi_approx 2 * area_half print(f半圓面積 {area_half:.16f}) print(f誤差估計(jì) {error_estimate:.2e}) import math print(fπ 近似值 {pi_approx:.16f}) print(fπ 真實(shí)值 {math.pi:.16f}) print(f絕對(duì)誤差 {abs(pi_approx - math.pi):.2e})我跑這段代碼得到的結(jié)果是半圓面積輸出1.5707963267948966π近似值3.141592653589793誤差在10的負(fù)16次方量級(jí)。這個(gè)精度已經(jīng)接近雙精度浮點(diǎn)數(shù)的極限了相當(dāng)驚人。重點(diǎn)看quad的第二個(gè)返回值error_estimate這是算法自己估算的絕對(duì)誤差。數(shù)值積分最妙的地方在于它不但給你結(jié)果還附帶一個(gè)可信度指標(biāo)。實(shí)際使用中我會(huì)習(xí)慣性把error_estimate打出來(lái)看一眼如果這個(gè)誤差估計(jì)很大說(shuō)明被積函數(shù)可能有問(wèn)題比如奇異點(diǎn)、震蕩劇烈或者區(qū)間選擇不當(dāng)。3.3 為什么quad這么準(zhǔn)自適應(yīng)算法揭秘也許你會(huì)想不就是求積分嗎我自己把區(qū)間切成10000份用矩形法累加也能算出差不多的π吧。沒(méi)錯(cuò)但效率差太多了。我做過(guò)一個(gè)簡(jiǎn)單測(cè)試用矩形法計(jì)算同樣的積分要達(dá)到誤差10的負(fù)8次方需要大約50萬(wàn)個(gè)月。而quad內(nèi)部采用自適應(yīng)Gauss-Kronrod積分它在函數(shù)變化緩的地方少采樣變化劇的地方多采樣。舉例來(lái)說(shuō)sqrt(1 - x2)在x接近±1時(shí)函數(shù)曲線幾乎垂直導(dǎo)數(shù)趨向無(wú)窮。普通均勻采樣在端點(diǎn)附近會(huì)很吃虧因?yàn)檫@一段函數(shù)變化太快均勻采樣的密度根本不夠。quad的聰明之處就在這里它會(huì)先在整體上采樣一組點(diǎn)然后檢查哪些區(qū)域的誤差貢獻(xiàn)大再對(duì)誤差大的區(qū)域遞歸加密采樣直到整個(gè)區(qū)間都滿足精度要求。這種“把力氣花在刀刃上”的策略正是自適應(yīng)算法的精髓。理解了這一點(diǎn)你就能解釋為什么同樣的積分quad耗時(shí)不到一毫秒精度卻碾壓手寫(xiě)矩形法幾個(gè)數(shù)量級(jí)。3.4 擴(kuò)展用quad做更復(fù)雜的圓面積計(jì)算既然標(biāo)題是“定積分求圓周率”我建議你做完基本版后順手做兩個(gè)擴(kuò)展練習(xí)。第一個(gè)擴(kuò)展是算任意半徑R的圓面積。只需要在half_circle里加入半徑參數(shù)用quad(half_circle, -R, R, args(R,))傳入。這里args的用法很實(shí)用后續(xù)做參數(shù)掃描時(shí)會(huì)經(jīng)常用到。def circle_upper(x, R): return np.sqrt(R ** 2 - x ** 2) for R in [1, 2, 0.5]: area, err quad(circle_upper, -R, R, args(R,)) print(f半徑{R}的圓面積: {2 * area:.10f}, 誤差估計(jì): {err:.1e})第二個(gè)擴(kuò)展是修改被積函數(shù)算半個(gè)橢圓的面積。橢圓方程x2/a2 y2/b2 1上半部分y b * sqrt(1 - x2/a2)積分區(qū)間[-a, a]。結(jié)果應(yīng)該是πab除以2。這個(gè)擴(kuò)展能讓你真正理解“定積分求面積”的本質(zhì)邊界函數(shù)求積分得到的是曲線下的面積如果把邊界函數(shù)改一下就能算各種稀奇古怪的圖形面積。4. 實(shí)戰(zhàn)中的坑與排查技巧4.1 FFT譜圖里那些奇怪的亮點(diǎn)和條紋第一次把頻譜用imshow顯示出來(lái)的時(shí)候很多人會(huì)嚇一跳為什么中心一個(gè)超亮的大斑點(diǎn)四周還有橫豎十字形的亮線中心大斑點(diǎn)其實(shí)是直流分量也就是全圖亮度的平均值。如果是一個(gè)全黑背景的圓盤(pán)圖直流分量會(huì)很高這是正常的。橫豎亮線通常來(lái)自圖像邊緣的不連續(xù)性圖像左右邊界、上下邊界在傅里葉眼里都是突變所以會(huì)在對(duì)應(yīng)的頻率軸上產(chǎn)生一條亮線。處理這類問(wèn)題一個(gè)經(jīng)典做法是給圖像做邊緣平滑或者加窗。不過(guò)在我們這個(gè)去噪項(xiàng)目里并不需要特別處理。只要你清楚這些現(xiàn)象的來(lái)源就不會(huì)被它們帶偏思路。4.2 逆變換結(jié)果出現(xiàn)復(fù)數(shù)怎么辦調(diào)用fft.ifft2之后返回結(jié)果其實(shí)是復(fù)數(shù)數(shù)組。這是因?yàn)楦↑c(diǎn)運(yùn)算誤差理論上實(shí)部才是真實(shí)圖像虛部會(huì)有極小的殘留。很多人忘了取實(shí)部直接plt.imshow復(fù)數(shù)數(shù)組結(jié)果要么報(bào)錯(cuò)要么顯示出一片亂碼。正確做法是調(diào)用.real取實(shí)部代碼里我已經(jīng)寫(xiě)了千萬(wàn)別省。還有一個(gè)常見(jiàn)坑是把ifftshift放在fftshift前面。正確的順序是fft2-fftshift- 濾波 -ifftshift-ifft2。如果你把ifftshift和fftshift搞混唯一的表現(xiàn)就是輸出圖像錯(cuò)亂這時(shí)候先檢查這一步。4.3 quad函數(shù)報(bào)錯(cuò)和異常情況處理quad最常見(jiàn)的報(bào)錯(cuò)是ValueError: Inf or NaN這通常是被積函數(shù)在積分區(qū)間內(nèi)出現(xiàn)無(wú)窮大。比如你把被積函數(shù)寫(xiě)成1 / sqrt(1 - x**2)在x ±1處分母為零函數(shù)值無(wú)窮大但恰恰這個(gè)積分是收斂的值等于π。quad其實(shí)能處理一部分端點(diǎn)奇異的情況但需要你顯式設(shè)置points參數(shù)把奇異點(diǎn)指出來(lái)。# 端點(diǎn)奇異但積分收斂的例子 integral, err quad(lambda x: 1 / np.sqrt(1 - x**2), -1, 1, points[-1, 1]) print(integral) # 接近 pi設(shè)置了points之后quad會(huì)在這些點(diǎn)附近加密采樣從而正確處理奇異行為。這個(gè)知識(shí)點(diǎn)在很多工程計(jì)算里都會(huì)用到值得記一下。另一個(gè)我遇到過(guò)的問(wèn)題是積分區(qū)間過(guò)大導(dǎo)致被積函數(shù)在中間出現(xiàn)劇烈震蕩。這時(shí)候quad可能報(bào)“Maximum number of subdivisions exceeded”的警告解決辦法是把大區(qū)間拆成多個(gè)小區(qū)間分別積分再把結(jié)果相加。實(shí)際工程里遇到這類問(wèn)題永遠(yuǎn)優(yōu)先考慮分區(qū)間而不是盲目調(diào)limit參數(shù)。4.4 常見(jiàn)問(wèn)題速查表癥狀可能原因解決辦法頻譜中心亮點(diǎn)巨大圖像直流分量過(guò)高正?，F(xiàn)象可做對(duì)數(shù)變換再顯示濾波后圖像有波紋使用了方形掩膜換成圓形掩膜降低振鈴逆變換輸出亂碼ifftshift位置錯(cuò)誤或沒(méi)取實(shí)部檢查順序加.real去噪后圖像過(guò)糊濾波半徑太小增大radius重新嘗試quad報(bào)Inf或NaN被積函數(shù)含奇異點(diǎn)用points參數(shù)標(biāo)示奇異點(diǎn)quad耗時(shí)長(zhǎng)且警告積分區(qū)間過(guò)大或震蕩劇烈手動(dòng)拆分積分區(qū)間分別計(jì)算這張表是我把實(shí)際踩過(guò)的坑整理出來(lái)的。做圖像處理時(shí)遇到視覺(jué)異常第一反應(yīng)先看頻譜和掩膜形狀對(duì)不對(duì)做數(shù)值積分時(shí)遇到報(bào)錯(cuò)第一反應(yīng)先把被積函數(shù)圖像畫(huà)出來(lái)。把函數(shù)可視化很多問(wèn)題一眼就能定位。我個(gè)人在這次實(shí)踐中最大的收獲是真正建立了“頻域思維”以前遇到圖像噪聲第一反應(yīng)是換濾波器窗口大小現(xiàn)在會(huì)先想到去頻譜里看一眼噪聲長(zhǎng)什么樣。這種思維方式的轉(zhuǎn)變遠(yuǎn)比記住幾個(gè)函數(shù)調(diào)用有價(jià)值。至于定積分求π看起來(lái)像個(gè)玩具但當(dāng)你意識(shí)到quad能處理帶奇異點(diǎn)的函數(shù)時(shí)你就會(huì)發(fā)現(xiàn)工程里很多復(fù)雜定積分都可以放心交給它前提是你要理解算法適應(yīng)的邊界在哪里。希望這兩個(gè)小項(xiàng)目能幫你打開(kāi)科學(xué)計(jì)算的大門(mén)。