周期函数的自身叠加后功率函数思考

周期函数的自身叠加后功率函数思考

根据傅里叶级数 任意周期为$\frac{2\pi }{w}$的函数f(t)
可以表示为 $f(t)=A_0+{\sum }_{n=1}^{\infty }{A}_{n}sin(nwt+{\psi }_n)$,n=1,2,3,…,$+\infty$

那么周期函数相移和经过幅值缩放后的函数可以表示为：
$f^{\prime }(t)=\lambda f(t+\tau )=B_0+{\sum }_{n=1}^{\infty }{B}_{n}sin(nwt+nw\tau +{\psi }_n)$,n=1,2,3,…,$+\infty$

周期函数相移和经过幅值缩放后与自身相加为 :$f(t)+\lambda f(t+\tau )$
根据信号与系统学，其功率函数为$\frac{w}{4\pi }{\int }_{-\frac{2\pi }{w}}^{\frac{2\pi }{w}}(f(t)+\lambda f(t+\tau ){)}^{2}dt$

先来研究一下$f(t)+\lambda f(t+\tau )$

它可以写成
$$\begin{gathered} f(t)+\lambda f(t+\tau )= \\ {A}_0+B_0+\sum _{n=1}^{\infty }\left(A_{n}sin(nwt+{\psi }_n)+B_{n}sin(nwt+nw\tau +{\psi }_n)\right) \end{gathered}$$

见下面的数学推导
$$\begin{array}{rl}Asin(nwt+\psi )+Bsin(nwt+nw\tau +\psi )& =\\ & =Asin(nwt+\psi )+Bsin(nwt+\psi )cos(nw\tau )+Bcos(nwt+\psi )sin(nw\tau )\\ & =(A+Bcos(nw\tau ))sin(nwt+\psi )+Bsin(nw\tau )cos(nwt+\psi )\\ & =\sqrt{A^2+B^2+2ABcos(nw\tau )}\left(\frac{A+Bcos(nw\tau )}{\sqrt{A^2+B^2+2ABcos(nw\tau )}}sin(nwt+\psi )+\frac{Bsin(nw\tau )}{\sqrt{A^2+B^2+2ABcos(nw\tau )}}cos(nwt+\psi )\right)\\ & =\sqrt{A^2+B^2+2ABcos(nw\tau )}\left(cos(\theta )sin(nwt+\psi )+sin(\theta )cos(nwt+\psi )\right)\\ & =\sqrt{A^2+B^2+2ABcos(nw\tau )}\ast sin(nwt+\psi +\theta )\end{array}$$
$A_{n}sin(nwt+{\psi }_n)+B_{n}sin(nwt+nw\tau +{\psi }_n)$又可以进行化简为：
$\sqrt{A_n^2+B_n^2+2A_n{B}_{n}cos(nw\tau )}\ast sin(nwt+{\psi }_n+{\theta }_n)$

所以
$$\begin{gathered} f(t)+\lambda f(t+\tau ) \\ =A_0+B_0+\sum _{n=1}^{\infty }\sqrt{A_n^2+B_n^2+2A_n{B}_{n}cos(nw\tau )}\ast sin(nwt+{\psi }_n+{\theta }_n) \\ =A_0^{\prime }+\sum _{n=1}^{\infty }{A}_n^{\prime }sin(nwt+{\psi }_n^{\prime }) \end{gathered}$$
其中 $A_0^{\prime }=A_0+B_0$
$A_n^{\prime }=\sqrt{A_n^2+B_n^2+2A_n{B}_{n}cos(nw\tau )}$
${\psi }_n^{\prime }={\psi }_n+{\theta }_n$

从上面的推导可知，$f(t)+\lambda f(t+\tau )$也是一个周期函数，其周期也是$\frac{2\pi }{w}$，
功率函数$\frac{w}{4\pi }{\int }_{-\frac{2\pi }{w}}^{\frac{2\pi }{w}}(f(t)+\lambda f(t+\tau ){)}^{2}dt$
相当于是从$A_0^{\prime },A_{n1}^{\prime }sin(n_{1}wt+{\psi }_{n1}^{\prime }),A_{n2}^{\prime }sin(n_{2}wt+{\psi }_{n2}^{\prime })...,A_n^{\prime }sin(nwt+{\psi }_n^{\prime })$任意抽出两个相乘后积分，

• 两个都是直流项
  $\frac{w}{4\pi }{\int }_{-\frac{2\pi }{w}}^{\frac{2\pi }{w}}{A}_0^{\prime }\ast A_0^{\prime }dt=\frac{w}{4\pi }{A}_0^{\prime }\ast A_0^{\prime }{|}_{-\frac{2\pi }{w}}^{\frac{2\pi }{w}}=A_0^{\prime }\ast A_0^{\prime }=(A_0+B_0{)}^2$

• 一个直流项 一个交流项
  $\frac{w}{4\pi }{\int }_{-\frac{2\pi }{w}}^{\frac{2\pi }{w}}{A}_0^{\prime }{A}_n^{\prime }sin(nwt+{\psi }_n^{\prime })dt=-\frac{w}{4\pi }{A}_0^{\prime }{A}_n^{\prime }\frac{1}{nw}cos(nwt+{\psi }_n^{\prime }){|}_{-\frac{2\pi }{w}}^{\frac{2\pi }{w}}=0$

• 两个都是交流项，且$n不等于m$
  $A_n^{\prime }sin(nwt+{\psi }_n^{\prime })A_m^{\prime }sin(mwt+{\psi }_m^{\prime })=\frac{A_n^{\prime }{A}_m^{\prime }}{2}\left(cos(nwt-mwt+{\psi }_n^{\prime }-{\psi }_m^{\prime })-cos(nwt+mwt+{\psi }_n^{\prime }+{\psi }_m^{\prime })\right)$
  由这个化简可得：
  $$\begin{array}{rl}\frac{w}{4\pi }{\int }_{-\frac{2\pi }{w}}^{\frac{2\pi }{w}}{A}_n^{\prime }sin(nwt+{\psi }_n^{\prime })A_m^{\prime }sin(mwt+{\psi }_m^{\prime })dt& =\\ & =\frac{w}{4\pi }\frac{A_n^{\prime }{A}_m^{\prime }}{2}\left(\frac{1}{nw-mw}sin(nwt-mwt+{\psi }_n^{\prime }-{\psi }_m^{\prime }){|}_{-\frac{2\pi }{w}}^{\frac{2\pi }{w}}-\frac{1}{nw+mw}sin(nwt+mwt+{\psi }_n^{\prime }+{\psi }_m^{\prime }){|}_{-\frac{2\pi }{w}}^{\frac{2\pi }{w}}\right)\\ & =0\end{array}$$

• 两个都是交流项，且$n等于m$
  $A_n^{\prime }sin(nwt+{\psi }_n^{\prime })A_m^{\prime }sin(mwt+{\psi }_m^{\prime })=\frac{A_n^{\prime }{A}_m^{\prime }}{2}\left(cos(nwt-mwt+{\psi }_n^{\prime }-{\psi }_m^{\prime })-cos(nwt+mwt+{\psi }_n^{\prime }+{\psi }_m^{\prime })\right)$
  =$\frac{A_n^{\prime }{A}_m^{\prime }}{2}\left(cos(0)-cos(2nwt+2{\psi }_n^{\prime })\right)$
  =$\frac{A_n^{\prime }{A}_n^{\prime }}{2}\left(1-cos(2nwt+2{\psi }_n^{\prime })\right)$

由这个化简可得：
$$\begin{array}{rl}\frac{w}{4\pi }{\int }_{-\frac{2\pi }{w}}^{\frac{2\pi }{w}}{A}_n^{\prime }sin(nwt+{\psi }_n^{\prime })A_m^{\prime }sin(mwt+{\psi }_m^{\prime })dt& =\\ & =\frac{w}{4\pi }\frac{A_n^{\prime }{A}_m^{\prime }}{2}\left(1{|}_{-\frac{2\pi }{w}}^{\frac{2\pi }{w}}-\frac{1}{2nw}sin(2nwt+2{\psi }_n^{\prime }){|}_{-\frac{2\pi }{w}}^{\frac{2\pi }{w}}\right)\\ & =\frac{w}{4\pi }\frac{A_n^{\prime }{A}_m^{\prime }}{2}\left(\frac{4\pi }{w}\right)\\ & =\frac{A_n^{\prime }{A}_n^{\prime }}{2}\\ & =\frac{A_n^2+B_n^2+2A_n{B}_{n}cos(nw\tau )}{2}\end{array}$$

综合四个可能，可知：
$$\begin{array}{rl}\frac{w}{4\pi }{\int }_{-\frac{2\pi }{w}}^{\frac{2\pi }{w}}(f(t)+\lambda f(t+\tau ){)}^{2}dt& =\\ & =(A_0+B_0{)}^2+\sum _{n=1}^{\infty }\frac{A_n^2+B_n^2+2A_n{B}_{n}cos(nw\tau )}{2}\\ & =(A_0+B_0{)}^2+\sum _{n=1}^{\infty }\frac{A_n^2+B_n^2}{2}+\sum _{n=1}^{\infty }{A}_n{B}_{n}cos(nw\tau )\end{array}$$

用ipython对上面推导出的公式进行验证:

import numpy as np
import math
from matplotlib import pyplot as plt 

%matplotlib inline

d1 = [0] * 1024
d2 = [0] * 1024
d3 = [0] * 1024
d4 = [0] * 1024
dsum = [0] * 1024 
p1 = [0] * 1024 
p2 = [0] * 1024 

psi = 0.9
a0 = 2.4
a1 = 2
a2 = 4
s = 0.5

for tau in range(1024):
    for i in range(1024):
        d1[i] = a1 * math.cos(3 * 2 * math.pi * i / 1024)
        d2[i] = a2 * math.cos(8 * 2 * math.pi * i / 1024 + psi)

        d3[i] = s * a1 * math.cos(3 * 2 * math.pi * (i + tau) / 1024)
        d4[i] = s * a2 * math.cos(8 * 2 * math.pi * (i + tau) / 1024 + psi)
        dsum[i] = a0 + d1[i] + d2[i] + s * a0 + d3[i] + d4[i]
        dsum[i] = dsum[i] ** 2

    p1[tau] = np.sum(dsum) / 1024
    p2[tau] = (a0 + s*a0)**2 + (a1**2 + a2**2 + (s*a1)**2 + (s*a2)**2)/2 + \
    a1*s*a1*math.cos(3 * 2 * math.pi * tau/1024) + a2*s*a2*math.cos(8 * 2 * math.pi * tau/1024)

    #print(p1,p2)

plt.figure(figsize=(15,10))
plt.plot(p1)
plt.plot(p2)

从输出的曲线来看，两段曲线完全重合

再来观察自叠加功率函数:
$$\begin{array}{rl}p(\tau )& =\frac{w}{4\pi }{\int }_{-\frac{2\pi }{w}}^{\frac{2\pi }{w}}(f(t)+\lambda f(t+\tau ){)}^{2}dt\\ & =(A_0+B_0{)}^2+\sum _{n=1}^{\infty }\frac{A_n^2+B_n^2}{2}+\sum _{n=1}^{\infty }{A}_n{B}_{n}cos(nw\tau )\end{array}$$
$cos(nw\tau )$总是在-1到1之间，
明显$\tau =0$时，$p(0)=(A_0+B_0{)}^2+{\sum }_{n=1}^{\infty }\frac{(A_n+B_n{)}^2}{2}$是$p(\tau )$的极大值。

---

---

光学相干长度相关的思考

参考了大神的文章
光的干涉原理
(其实主要是为了搞懂干涉显微的大概原理 干涉术：白光扫描干涉显微镜 )
里面推导有点复杂，但是大概的要推导出相干长度的计算公式，需要假设光源的光谱为一条幅值为$I_0$的连续直线，相干条纹是由这样的两个相同光源叠加生成的。

我们设光的三角函数波动模型为: $f(t)={\sum }_{n_a}^{n_b}{i}_{0}sin(nwt+{\psi }_n)$,(没有直流项的傅里叶级数，$w$趋向于无穷小,n_b大于n_a)
若C为光速，波长为$\lambda$，$nw=2\pi \frac{C}{\lambda }$，光谱起始波长为${\lambda }_b=2\pi \frac{C}{n_{b}w}$，光谱结束波长为${\lambda }_a=2\pi \frac{C}{n_{a}w}$。
而两个相同光源叠加的功率函数可表示为干涉条纹函数，根据上面功率函数相关推导有：
$$\begin{array}{rl}p(\tau )& =\frac{w}{4\pi }{\int }_{-\frac{2\pi }{w}}^{\frac{2\pi }{w}}(f(t)+f(t+\tau ){)}^{2}dt\\ & =\sum _{n_a}^{n_b}\frac{i_0^2+i_0^2}{2}+\sum _{n_a}^{n_b}{i}_0^{2}cos(nw\tau )\\ & =\sum _{n_a}^{n_b}{i}_0^2+\sum _{n_a}^{n_b}{i}_0^{2}cos(nw\tau )\end{array}$$
由于$i_0$是一个常数，而$w$比较小，所以常数和cos求和会约等于积分：
$$\begin{array}{rl}p(\tau )& =\sum _{n_a}^{n_b}{i}_0^2+\sum _{n_a}^{n_b}{i}_0^{2}cos(nw\tau )\\ & \approx {\int }_{n_a}^{n_b}{i}_0^{2}dn+{\int }_{n_a}^{n_b}{i}_0^{2}cos(nw\tau )dn\\ & \approx (n_b-na)i_0^2+\frac{i_0^{2}sin(nw\tau )}{w\tau }{|}_{n_a}^{n_b}\\ & \approx (n_b-na)i_0^2+\frac{i_0^{2}sin(n_{b}w\tau )-i_0^{2}sin(n_{a}w\tau )}{w\tau }\\ & \approx (n_b-na)i_0^2(1+\frac{sin(n_{b}w\tau )-sin(n_{a}w\tau )}{(n_b-na)w\tau })\end{array}$$

假设光源的中心波长为${\lambda }_c$，则有对应的中心谱数$n_c=\frac{n_a+n_b}{2}$，${\lambda }_c=2\pi \frac{C}{n_{c}w}$

假设光源的光谱带度为$\Delta \lambda ={\lambda }_a-{\lambda }_b$，则有对应的谱数范围$\Delta n=n_b-n_a$,
$n_{a}w=(n_c-\Delta n/2)w$
$n_{b}w=(n_c+\Delta n/2)w$
$\Delta \lambda =2\pi (\frac{C}{n_{a}w}-\frac{C}{n_{b}w})$
所以：：
$$\begin{array}{rl}p(\tau )& \approx (n_b-na)i_0^2(1+\frac{sin(n_{b}w\tau )-sin(n_{a}w\tau )}{(n_b-na)w\tau })\\ & \approx \Delta n{i}_0^2(1+\frac{sin(n_{c}w\tau +\Delta nw\tau /2)-sin(n_{c}w\tau -\Delta nw\tau /2)}{\Delta nw\tau })\\ & \approx \Delta n{i}_0^2(1+\frac{2cos(n_{c}w\tau )sin(\Delta nw\tau /2)}{\Delta nw\tau })\\ & \approx \Delta n{i}_0^2(1+\frac{cos(n_{c}w\tau )sin(\Delta nw\tau /2)}{\Delta nw\tau /2})\end{array}$$

用ipython进行上面公式的验证:

spectralf = [0+0j]*1024
cowave = [0]*1024
P_wave = [0]*1024

s = 100
e = 110
center = (s+e)/2
delta = e-s
omiga = math.pi / (e - s + 1)
for i in range(1024):
    if s<=i<=e:
        spectralf[i] = 512+0j  #乘512是由于IFFT是以N*X(k)来作为频谱幅值来算的
    if  1024 - e <=i<= 1024 - s:
        spectralf[i] = 512+0j


# spectralf = [0+0j]*1024
# spectralf[1] = 1 + 0j
# spectralf[1023] = 1 + 0j
wwave = np.fft.ifft(spectralf)
realwwave = np.real(wwave)
for i in range(1024):
    sumd = 0
    for j in range(1024):
        t = j + i
        if t>1023:
            t=t%1024
        sumd = sumd + (realwwave[j] + realwwave[t])**2
    cowave[i] = (2*math.pi/1024/2/math.pi)* sumd #这里积分只用了一个周期(0~2*pi/w) 所以要除以(2*pi/w - 0)
    
    tau = i % 512
    if i>511:
        tau = 512-tau #强制相移不超过半周期，超过相当于负相移
    if tau>0:    
        P_wave[i]=delta * (1 + \
                           math.cos(center * 2 * math.pi * tau/1024)* \
                           math.sin(delta * 2 * math.pi * tau/1024/2)/  \
                           (delta * 2 * math.pi * tau/1024/2) \
                          )

P_wave[0] = 1

plt.figure(figsize=(15,10))
plt.plot(cowave)
plt.plot(P_wave)
i = int(2*1024/(delta * 2 ) )
plt.plot([1*delta*2]*(i)+[0])

上图中 蓝色波形为原始波形直接运算叠加的干涉强度(功率函数)波形图，橙色波形是用推导的公式运算干涉强度，可以看到，橙色与蓝色十分接近，而且计算出的第一个包络过零点与蓝色的波形一致。

---

---

知乎光的干涉原理 上说干涉条纹可见度为:
$$V=\left|\frac{sin(\Delta nw\tau /2)}{\Delta nw\tau /2}\right|$$
但我搞不懂这是从那里来的-_-!，从其它的方向来分析吧。
$\Delta n{i}_0^2(1+\frac{cos(n_{c}w\tau )sin(\Delta nw\tau /2)}{\Delta nw\tau /2})$这个式子其实可以看到

其振动部分主要是$\frac{cos(n_{c}w\tau )sin(\Delta nw\tau /2)}{\Delta nw\tau /2}$
而$n_c=\frac{n_a+n_b}{2}$ 和 $\frac{\Delta n}{2}=\frac{n_b-n_a}{2}$
明显有$\frac{\Delta n}{2}/n_c=\frac{n_b-n_a}{n_a+n_b}<1$ ,即 $n_c>\frac{\Delta n}{2}$

再来看看$\frac{sinx}{x}$的特性:
x趋向于0时，sinx的一阶导数为cosx,也就是要说x比较小时，
$sinx\approx sin(0)+cos(0)(x-0)=x$
而当x比0大得多时，相当于是振动波形sinx乘以了一个$\frac{1}{x}$的包络，
以ipython 演示$\frac{sinx}{x}$波形

wave = [0]*1024
wave11 = [0]*1024
wave11[0] = 1
wave[0] = 1
for i in range(1,1024):
    tt =4 * 2 * math.pi * i / 1024
    wave[i] = math.sin(tt)/tt
    wave11[i] = math.cos(tt*18.51)*wave[i]

plt.figure(figsize=(15,10))
plt.plot([0]*1024)
plt.plot(wave)
plt.figure(figsize=(15,10))
plt.plot([0]*1024)
plt.plot(wave11)
i = int(1024/(4 * 2 ) )
plt.plot([1]*i+[0])

振动部分$\frac{cos(n_{c}w\tau )sin(\Delta nw\tau /2)}{\Delta nw\tau /2}$相当于是在上图的波形乘以更快速度振动$cos(n_{c}w\tau )$项，也就是下面的波形:

从图中就可以看到，振动部分$\frac{cos(n_{c}w\tau )sin(\Delta nw\tau /2)}{\Delta nw\tau /2}$第一个比较明显的包络过零点应该就是最大的相干相移${\tau }_{max}$了，也就是$sin(\Delta nw{\tau }_{max}/2)=0,\Delta nw{\tau }_{max}/2=\pi$
对${\tau }_{max}$进行变换
$$\begin{array}{rl}{\tau }_{max}& =\frac{2\pi }{\Delta nw}\\ & =\frac{2\pi }{\Delta nw}\\ & =\frac{2\pi }{(n_b-n_a)w},相当于\frac{T}{(n_b-n_a)}\end{array}$$
相干长度$L=C{\tau }_{max}$，(光速度乘以最大相干相移时长
波长公式:$\lambda =2\pi \frac{C}{nw}$
所以：：
$$\begin{array}{rl}L& =C{\tau }_{max}\\ & =C\frac{2\pi }{(n_b-n_a)w}\\ & =\frac{2\pi C}{n_{b}w-n_{a}w}\\ & =\frac{2\pi C}{2\pi (\frac{C}{{\lambda }_b}-\frac{C}{{\lambda }_a})}\\ & =\frac{1}{\frac{1}{{\lambda }_b}-\frac{1}{{\lambda }_a}}\\ & =\frac{1}{\frac{1}{{\lambda }_b}-\frac{1}{{\lambda }_a}}\\ & =\frac{{\lambda }_a{\lambda }_b}{{\lambda }_a-{\lambda }_b}\\ & =\frac{{\lambda }_a{\lambda }_b}{\Delta \lambda }\end{array}$$
又由于现实中干涉测量用的光源光谱起始波长与结束波长与中心波长${\lambda }_c$相差比较小
所以
$$\begin{array}{rl}L& =\frac{{\lambda }_a{\lambda }_b}{\Delta \lambda }\\ & \approx \frac{{\lambda }_c^2}{\Delta \lambda }\end{array}$$