UA OPTI512R 傅立叶光学导论9 卷积基础

卷积的图示
卷积的解析计算

之前介绍LSI时，我们提到过如果输入为δ(x)\delta(x)δ(x)时，输出为脉冲响应函数h(x)=L[δ(x)]h(x)=\mathcal{L}[\delta(x)]h(x)=L[δ(x)]则输出为f(x)f(x)f(x)时，输出为输入与脉冲响应函数的卷积
g(x)=∫−∞+∞f(y)h(x−y)dy=(f∗h)(x)g(x)=\int_{-\infty}^{+\infty} f(y)h(x-y)dy=(f*h)(x)g(x)=∫−∞+∞f(y)h(x−y)dy=(f∗h)(x)

做一个简单的自变量变换z=x−yz=x-yz=x−y，则
∫−∞+∞f(y)h(x−y)dy=∫−∞+∞h(z)f(x−z)dz=(h∗f)(x)\int_{-\infty}^{+\infty} f(y)h(x-y)dy=\int_{-\infty}^{+\infty}h(z)f(x-z)dz = (h*f )(x)∫−∞+∞f(y)h(x−y)dy=∫−∞+∞h(z)f(x−z)dz=(h∗f)(x)

所以f∗hf*hf∗h与h∗fh*fh∗f这两种记法完全一样。这一讲介绍一些卷积的基础知识。

卷积的图示

假设
f(x)=rect(x−12)h(x)=ramp(x)rect(x−1.53)f(x) = rect \left( \frac{x-1}{2} \right) \\ h(x) = ramp(x)rect \left( \frac{x-1.5}{3} \right)f(x)=rect(2x−1)h(x)=ramp(x)rect(3x−1.5)

我们尝试图示
(f∗h)(x)=∫−∞+∞f(y)h(x−y)dy(f*h)(x)=\int_{-\infty}^{+\infty} f(y)h(x-y)dy(f∗h)(x)=∫−∞+∞f(y)h(x−y)dy

计算重叠区域面积的过程图示如下

卷积的解析计算

现在尝试直接计算
(f∗h)(x)=∫−∞+∞f(y)h(x−y)dy=∫−∞+∞rect(y−12)ramp(x−y)rect((x−y)−1.53)dy={0,x<0∫0x(x−y)dy=x22,0≤x<2∫02(x−y)dy=2(x−1),2≤x<3∫x−32(x−y)dy=−x2+4x+52,3≤x<50,x≥5\begin{aligned} (f*h)(x)& =\int_{-\infty}^{+\infty} f(y)h(x-y)dy \\ & = \int_{-\infty}^{+\infty}rect\left( \frac{y-1}{2} \right)ramp(x-y)rect \left( \frac{(x-y)-1.5}{3} \right)dy \\ & = \begin{cases} 0 ,x < 0 \\ \int_0^x (x-y)dy = \frac{x^2}{2}, 0 \le x < 2 \\ \int_0^2 (x-y)dy = 2(x-1),2 \le x < 3 \\ \int_{x-3}^2 (x-y)dy = \frac{-x^2+4x+5}{2},3 \le x < 5 \\ 0, x \ge 5 \end{cases}\end{aligned}(f∗h)(x)=∫−∞+∞f(y)h(x−y)dy=∫−∞+∞rect(2y−1)ramp(x−y)rect(3(x−y)−1.5)dy=⎩⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎧0,x<0∫0x(x−y)dy=2x2,0≤x<2∫02(x−y)dy=2(x−1),2≤x<3∫x−32(x−y)dy=2−x2+4x+5,3≤x<50,x≥5

怎么说呢，要解析地算出卷积的表达式其实就是算积分或者求和/求级数，如果impulse response和input的形式都很简单，解析计算还可以，但函数形式比较复杂的话解析计算就行不通了，后续在介绍完卷积的性质后会介绍计算卷积的数值方法。

例1 计算下面两个函数的卷积
f(x)=rect(x−12),h(x)=2rect(x−0.5)f(x)=rect\left( \frac{x-1}{2} \right),h(x)=2rect(x-0.5)f(x)=rect(2x−1),h(x)=2rect(x−0.5)

(f∗h)(x)=∫−∞+∞f(y)h(x−y)dy=2∫−∞+∞rect(y−12)rect((x−y)−0.5)dy=2∫−∞+∞χ(0,2)(y)χ(x−1,x)(y)dy=2∫−∞+∞χ(0,2)∩(x−1,x)(y)dy={0,x<02∫0xdy=2x,0≤x<12∫x−1xdy=2,1≤x<22∫x−12dy=6−2x,2≤x<30,x≥3\begin{aligned} (f*h)(x) & = \int_{-\infty}^{+\infty} f(y) h(x-y)dy \\ & = 2\int_{-\infty}^{+\infty} rect\left( \frac{y-1}{2} \right) rect((x-y)-0.5)dy \\ & = 2\int_{-\infty}^{+\infty} \chi_{(0,2)}(y) \chi_{(x-1,x)}(y)dy \\ & =2\int_{-\infty}^{+\infty} \chi_{(0,2) \cap (x-1,x)}(y)dy \\ & = \begin{cases} 0, x < 0 \\ 2 \int_0^x dy = 2x, 0 \le x < 1 \\ 2\int_{x-1}^x dy = 2, 1 \le x < 2 \\ 2 \int_{x-1}^2 dy = 6-2x, 2 \le x < 3 \\ 0 , x \ge 3 \end{cases}\end{aligned}(f∗h)(x)=∫−∞+∞f(y)h(x−y)dy=2∫−∞+∞rect(2y−1)rect((x−y)−0.5)dy=2∫−∞+∞χ(0,2)(y)χ(x−1,x)(y)dy=2∫−∞+∞χ(0,2)∩(x−1,x)(y)dy=⎩⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎧0,x<02∫0xdy=2x,0≤x<12∫x−1xdy=2,1≤x<22∫x−12dy=6−2x,2≤x<30,x≥3

例2 卷积与Fresnel衍射

考虑二元函数作为LSI的输入：
g(x1,x2)=L[f(x1,x2)]g(x_1,x_2) = \mathcal{L}[f(x_1,x_2)]g(x1,x2)=L[f(x1,x2)]

假设δ(x1,x2)\delta(x_1,x_2)δ(x1,x2)作为输入时脉冲响应函数为h(x1,x2)h(x_1,x_2)h(x1,x2)，则
g(x1,x2)=∬−∞+∞f(y1,y2)h(x1−y1,x2−y2)dy1dy2=(f∗h)(x1,x2)g(x_1,x_2) = \iint_{-\infty}^{+\infty} f(y_1,y_2) h(x_1-y_1,x_2-y_2)dy_1dy_2 =( f*h)(x_1,x_2)g(x1,x2)=∬−∞+∞f(y1,y2)h(x1−y1,x2−y2)dy1dy2=(f∗h)(x1,x2)

这类结果可以自然推广到nnn元函数。考虑x=(x1,x2,⋯,xn)′\textbf x=(x_1,x_2,\cdots,x_n)'x=(x1,x2,⋯,xn)′,
g(x)=L[f(x)]g(\textbf x) = \mathcal{L}[f(\textbf x)]g(x)=L[f(x)]

假设δ(x)\delta(\textbf x)δ(x)作为输入时脉冲响应函数为h(x)h(\textbf x)h(x)，则
g(x)=∫−∞+∞f(y)h(x−y)dy=(f∗h)(x)g(\textbf x) = \int_{-\infty}^{+\infty} f(\textbf y) h(\textbf x-\textbf y)d\textbf y =( f*h)(\textbf x)g(x)=∫−∞+∞f(y)h(x−y)dy=(f∗h)(x)

考虑波前为ui(x,y)u_i(x,y)ui(x,y)的光，假设材料透光性用P(x,y)P(x,y)P(x,y)描述，比如如果是圆孔的话可以用cylinder function，经过材料后，波前变为u1(x,y)=ui(x,y)P(x,y)u_1(x,y)=u_i(x,y)P(x,y)u1(x,y)=ui(x,y)P(x,y)

根据Rayleigh-Sommerfeld diffraction integral的Kirchoff近似，将波前u1u_1u1上的每一点(x′,y′,0)(x',y',0)(x′,y′,0)视为新的光源，传播到(x,y,z)(x,y,z)(x,y,z)时，波前为
u2(x,y)=∬−∞+∞u1(x′,y′)zjλrejkrrdx′dy′u_2(x,y)=\iint_{-\infty}^{+\infty}u_1(x',y') \frac{z}{j \lambda r} \frac{e^{jkr}}{r}dx'dy'u2(x,y)=∬−∞+∞u1(x′,y′)jλrzrejkrdx′dy′

使用Fresnel衍射的近似方法（Fresnel衍射的含义是我们观察u2u_2u2的地方与透光材料之间的距离相比透光材料的孔径非常大，所以在z/(jλr)z/(j \lambda r)z/(jλr)中，取z≈rz \approx rz≈r，然后用前两阶Taylor展开近似ejkrr\frac{e^{jkr}}{r}rejkr），
u2(x,y)=ejkzjλz∬P(x,y)ui(x′,y′)ejπλz[(x−x′)2+(y−y′)2]dx′dy′=(u1∗h)(x,y),h=ejπλz(x2+y2)\begin{aligned}u_2(x,y) & =\frac{e^{jkz}}{j \lambda z} \iint_{P(x,y)}u_i(x',y')e^{j \frac{\pi}{\lambda z}[(x-x')^2+(y-y')^2]}dx'dy' \\ & = (u_1*h)(x,y),h=e^{j \frac{\pi}{\lambda z}(x^2+y^2)}\end{aligned}u2(x,y)=jλzejkz∬P(x,y)ui(x′,y′)ejλzπ[(x−x′)2+(y−y′)2]dx′dy′=(u1∗h)(x,y),h=ejλzπ(x2+y2)

也就是说，当我们并不关注这个衍射系统的设计，只想知道输入的波形与输出的波形之间的关系时，可以把它看成一个LSI，其中hhh被称为Fresnel Diffraction Impulse Response。