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