时域OCT(TD-OCT) 我的理解

时域OCT(TD-OCT) 我的理解

由于本人是先接触SS-OCT的 所以对TD-OCT有点不理解,百度查找了一些资料总结一下。

TD-OCT 光学系统跟SS-OCT差不多,如下图
在这里插入图片描述
引用知乎上的原理解释:https://zhuanlan.zhihu.com/p/46452072
在这里插入图片描述时域OCT[1]SLD发出的低相干宽谱光在单模光纤中传递(如果使用多模光纤,模式色散会展宽宽光谱光,导致OCT的轴向分辨率下降),经过一个50:50的光纤耦合器(coupler)后,各有50%的光分别进入到测量臂和参考臂中照射到样品和参考镜上,来自样品的背向散射光和来自参考镜的反射光在耦合器上合并后被探测器(detector)收集。需要注意的是,只有当两臂的光程差在宽光谱光源相干长度的范围内(即光程差﹤光源的相干长度),两束光才可以发生干涉,这也使得OCT可以分辨出样品背向散射光的轴向长度分布。因此可以在参考臂利用步进电机在轴向平移参考镜(Reference),实现轴向扫描,得到样品的轴向深度分布图。同时在样品臂对光束进行横向扫描,可以得到样品的二维横截面图像,进而产生三维立体图像。
上述OCT技术是基于时间域的低相干干涉测量技术,通常称之为“时域OCT(Time Domain OCT, TD-OCT)”。时域OCT技术出现之后,由于其可以兼顾较高的成像深度和分辨率,迅速得到了研究者的广泛关注。但是在之后的研究中发现时域OCT在进行轴向扫描时,由于步进电机的限制,会致测量速度受到限制,机械装置的引入也会使系统的稳定性大打折扣。针对这些缺点,傅里叶域OCT(Fourier Domain OCT, FD-OCT)应运而生。基于测量臂与参考臂的干涉光谱与样品不同深度的背向散射光强度信息恰好是一对傅里叶变换对的关系这一原理,傅里叶域OCT避免了在测量时需要改变参考臂光程长度的问题,而是将代表着样品不同深度信息的光束统一收集后,通过傅里叶(逆)变换即可解算出样品深度的强度信息,实现了样品轴向信息的并行获取,从而省却了轴向扫描。

这里的"SLD发出的低相干宽谱光" 我认为是关键信息
找了一圈百度也没有找到SLD(超辐射发光二极管),能发出什么波形的光(不正确的说法),直到看到一段话:
如果光源在频域是宽带,那么在空间域就是窄带
这段话也是出于知乎的 https://zhuanlan.zhihu.com/p/159419026
然后又看到了大佬的文章 http://www.elecfans.com/d/845686.html里面的两个插图

在这里插入图片描述在这里插入图片描述
我就在这文章里就武断地认为SLD光源发出的光的波的特性如下:

  • 在时间上是分段的波形
  • 形状上估计差不多是中心频率正弦波乘以正态分布的一个波形吧
  • 如果对波形进行傳里叶变换,变换的结果相当于在中心频率上有一个差不多正态分布形状的频谱

这些特性是我个人猜测的,并没有查找到实际根据。

根据上面找到的信息,我自己用ipython做一下原理上的模拟如下

生成模拟SLD光源的光波形

import numpy as np
import math
import cmath

%pylab inline
x = [0] * 512

o = 10 #标准差越小 波形越窄 对应的频谱越宽
u = 255
centerfreq = 80
for i in range(512):
    t = 1 / math.sqrt(2 * math.pi ) * o
    t1 = math.exp(-(i - u)**2 / (2*o**2))
    x[i] = t*t1*math.sin(2 * math.pi * centerfreq * i / 511)

plt.plot(x)

#x = np.random.rand(64)
xfft = np.fft.fft(x)
xfft = np.abs(xfft)

plt.plot(xfft[:len(xfft)//2])
print("傳里叶变换后的中心频率 =",xfft.argmax())

在这里插入图片描述
傳里叶变换后的中心频率 = 80
在这里插入图片描述
也可以从正态分布的频域生成波形

print("从宽带谱域反过来生成波形")
specwave = [0]*512
specwidth = 10 #标准差越小 波形越窄 对应的频谱越宽
centerfreq = 80
for i in range(256):
    t = 1 / math.sqrt(2 * math.pi ) * specwidth
    t1 = math.exp(-(i - centerfreq)**2 / (2*specwidth**2))
    specwave[i] = t*t1
    if i>0 :  #频域理论上应该是对称的 (幅值相等 相位相反)
        specwave[512-i] = specwave[i]
specwave[256] = 0   
ifff = np.fft.ifft(specwave)
ifff = np.concatenate((ifff[-256:],ifff[:-256]))
plt.figure()
plt.text(0,0.2,"Spectral")
plt.plot(specwave)
plt.figure()
plt.text(0,0.2,"wave")
plt.plot(ifff.real)

从宽带谱域反过来生成波形
在这里插入图片描述在这里插入图片描述

从一个SLD光源 获取的波长光谱 生成波形


specdata_wlen = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
                 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 
                 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
                 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 
                 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 
                 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 
                 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 
                 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 
                 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
                 0.0, 0.0, 0.0, 0.0, 0.0018315018315018315, 0.0018315018315018315, 
                 0.0018315018315018315, 0.0018315018315018315, 0.0018315018315018315,
                 0.0018315018315018315, 0.0018315018315018315, 0.0, 0.0, 0.0, 0.0, 
                 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.003663003663003663, 
                 0.003663003663003663, 0.003663003663003663, 0.003663003663003663,
                 0.003663003663003663, 0.003663003663003663, 0.003663003663003663, 
                 0.003663003663003663, 0.003663003663003663, 0.003663003663003663, 0.003663003663003663, 
                 0.007326007326007326, 0.007326007326007326,
                 0.007326007326007326, 0.007326007326007326, 0.007326007326007326, 
                 0.007326007326007326, 0.007326007326007326, 0.007326007326007326, 0.01098901098901099,
                 0.01098901098901099, 0.01098901098901099, 0.01098901098901099, 0.01098901098901099, 
                 0.01282051282051282, 0.01282051282051282, 0.01282051282051282, 0.016483516483516484,
                 0.016483516483516484, 0.016483516483516484, 0.016483516483516484, 0.020146520146520148, 
                 0.020146520146520148, 0.020146520146520148, 0.023809523809523808, 0.023809523809523808, 0.027472527472527472, 
                 0.027472527472527472, 0.027472527472527472, 0.027472527472527472, 0.029304029304029304, 0.031135531135531136, 
                 0.03296703296703297, 0.03663003663003663, 0.03663003663003663, 0.040293040293040296, 0.040293040293040296, 
                 0.04395604395604396, 0.04395604395604396, 0.045787545787545784, 0.04945054945054945, 0.04945054945054945, 
                 0.05311355311355311, 0.05311355311355311, 0.056776556776556776, 0.056776556776556776, 0.06227106227106227,
                 0.0641025641025641, 0.06593406593406594, 0.07326007326007326, 0.07326007326007326, 0.07509157509157509, 
                 0.07509157509157509, 0.08241758241758242, 0.08241758241758242, 0.08791208791208792, 0.09523809523809523, 
                 0.09523809523809523, 0.10256410256410256, 0.10256410256410256, 0.10805860805860806, 0.10805860805860806,
                 0.11538461538461539, 0.11538461538461539, 0.12087912087912088, 0.1282051282051282, 0.1282051282051282, 
                 0.13736263736263737, 0.13736263736263737, 0.1446886446886447, 0.1446886446886447, 0.15384615384615385, 
                 0.15384615384615385, 0.163003663003663, 0.17399267399267399, 0.17399267399267399, 0.18315018315018314, 
                 0.18315018315018314, 0.19230769230769232, 0.19413919413919414, 0.2032967032967033, 0.21245421245421245, 
                 0.21245421245421245, 0.22344322344322345, 0.22344322344322345, 0.2326007326007326, 0.2326007326007326, 
                 0.24175824175824176, 0.24358974358974358, 0.25457875457875456, 0.26556776556776557, 0.26556776556776557, 
                 0.2783882783882784, 0.2783882783882784, 0.2875457875457875, 0.2875457875457875, 0.29853479853479853, 
                 0.3076923076923077, 0.30952380952380953, 0.32051282051282054, 0.32051282051282054, 0.32967032967032966, 
                 0.32967032967032966, 0.34065934065934067, 0.34065934065934067, 0.3534798534798535, 0.3626373626373626, 
                 0.3626373626373626, 0.3717948717948718, 0.37362637362637363, 0.38278388278388276, 0.38278388278388276, 
                 0.39194139194139194, 0.39377289377289376, 0.40293040293040294, 0.41208791208791207, 0.41208791208791207, 
                 0.42124542124542125, 0.42124542124542125, 0.42857142857142855, 0.42857142857142855, 0.43772893772893773, 
                 0.44505494505494503, 0.44505494505494503, 0.45054945054945056, 0.45054945054945056, 0.45787545787545786, 
                 0.45787545787545786, 0.4633699633699634, 0.4652014652014652, 0.4706959706959707, 0.47619047619047616, 
                 0.47802197802197804, 0.47985347985347987, 0.47985347985347987, 0.48717948717948717, 0.48717948717948717,
                 0.4908424908424908, 0.4926739926739927, 0.4926739926739927, 0.49633699633699635, 0.49633699633699635, 0.5, 
                 0.5, 0.5, 0.5, 0.5036630036630036, 0.5036630036630036, 0.5036630036630036, 0.5036630036630036, 0.5036630036630036,
                 0.5036630036630036, 0.5036630036630036, 0.5036630036630036, 0.5036630036630036, 0.5036630036630036, 0.5036630036630036, 
                 0.5018315018315018, 0.5, 0.5, 0.5, 0.49633699633699635, 0.49633699633699635, 0.49633699633699635, 0.49633699633699635, 
                 0.4926739926739927, 0.4908424908424908, 0.4908424908424908, 0.48717948717948717, 0.48717948717948717, 0.48717948717948717,
                 0.48717948717948717, 0.4835164835164835, 0.4816849816849817, 0.47985347985347987, 0.47802197802197804, 0.47802197802197804, 
                 0.47435897435897434, 0.47435897435897434, 0.47435897435897434, 0.47435897435897434, 0.4706959706959707, 0.46703296703296704,
                 0.46703296703296704, 0.4633699633699634, 0.4633699633699634, 0.46153846153846156, 0.46153846153846156, 0.45787545787545786, 
                 0.45787545787545786, 0.4542124542124542, 0.4523809523809524, 0.45054945054945056, 0.44871794871794873, 0.4468864468864469,
                 0.44505494505494503, 0.44505494505494503, 0.4413919413919414, 0.4413919413919414, 0.4413919413919414, 0.43772893772893773,
                 0.43772893772893773, 0.4340659340659341, 0.4340659340659341, 0.43223443223443225, 0.43223443223443225, 0.43223443223443225,
                 0.42857142857142855, 0.42857142857142855, 0.4249084249084249, 0.4249084249084249, 0.4249084249084249, 0.4249084249084249,
                 0.42124542124542125, 0.4194139194139194, 0.4175824175824176, 0.4175824175824176, 0.4175824175824176, 0.4157509157509158, 
                 0.4157509157509158, 0.4157509157509158, 0.4157509157509158, 0.41208791208791207, 0.41208791208791207, 0.41208791208791207,
                 0.4084249084249084, 0.4084249084249084, 0.4084249084249084, 0.4084249084249084, 0.4084249084249084, 0.4084249084249084,
                 0.40476190476190477, 0.40476190476190477, 0.40476190476190477, 0.40476190476190477, 0.40476190476190477, 0.40476190476190477,
                 0.40476190476190477, 0.40476190476190477, 0.40476190476190477, 0.40476190476190477, 0.40476190476190477, 0.40476190476190477,
                 0.40476190476190477, 0.40476190476190477, 0.40476190476190477, 0.40476190476190477, 0.40476190476190477, 0.40476190476190477,
                 0.40476190476190477, 0.40476190476190477, 0.40476190476190477, 0.40476190476190477, 0.40476190476190477, 0.40476190476190477,
                 0.40476190476190477, 0.40476190476190477, 0.40476190476190477, 0.40476190476190477, 0.4065934065934066, 0.4084249084249084, 
                 0.4084249084249084, 0.4084249084249084, 0.4084249084249084, 0.4084249084249084, 0.4084249084249084, 0.4084249084249084, 
                 0.41208791208791207, 0.41208791208791207, 0.41208791208791207, 0.41208791208791207, 0.41208791208791207, 0.4157509157509158,
                 0.4157509157509158, 0.4157509157509158, 0.4157509157509158, 0.4175824175824176, 0.4175824175824176, 0.4175824175824176,
                 0.42124542124542125, 0.42124542124542125, 0.4249084249084249, 0.4249084249084249, 0.4249084249084249, 0.4267399267399267, 
                 0.42857142857142855, 0.42857142857142855, 0.42857142857142855, 0.43223443223443225, 0.43223443223443225, 0.4340659340659341,
                 0.4340659340659341, 0.4340659340659341, 0.43772893772893773, 0.43772893772893773, 0.4413919413919414, 0.4413919413919414, 
                 0.44505494505494503, 0.44505494505494503, 0.44505494505494503, 0.44505494505494503, 0.4468864468864469, 0.45054945054945056,
                 0.45054945054945056, 0.4542124542124542, 0.4542124542124542, 0.45787545787545786, 0.45787545787545786, 0.46153846153846156,
                 0.46153846153846156, 0.4633699633699634, 0.46703296703296704, 0.46703296703296704, 0.46703296703296704, 0.46703296703296704,
                 0.4706959706959707, 0.4706959706959707, 0.47435897435897434, 0.47435897435897434, 0.47619047619047616, 0.47985347985347987, 
                 0.47985347985347987, 0.4835164835164835, 0.4835164835164835, 0.4908424908424908, 0.4908424908424908, 0.4926739926739927, 
                 0.49633699633699635, 0.49633699633699635, 0.5, 0.5, 0.5036630036630036, 0.5036630036630036, 0.5073260073260073,
                 0.5091575091575091, 0.5091575091575091, 0.5128205128205128, 0.5128205128205128, 0.5164835164835165, 0.5164835164835165,
                 0.5201465201465202, 0.5201465201465202, 0.521978021978022, 0.5256410256410257, 0.5256410256410257, 0.532967032967033, 
                 0.532967032967033, 0.5366300366300366, 0.5366300366300366, 0.5384615384615384, 0.5421245421245421, 0.5421245421245421, 
                 0.5494505494505495, 0.5494505494505495, 0.5512820512820513, 0.5531135531135531, 0.554945054945055, 0.5567765567765568, 
                 0.5586080586080586, 0.5622710622710623, 0.5622710622710623, 0.5677655677655677, 0.5695970695970696, 0.5714285714285714, 
                 0.5714285714285714, 0.575091575091575, 0.5769230769230769, 0.5787545787545788, 0.5842490842490843, 0.5842490842490843,
                 0.5879120879120879, 0.5879120879120879, 0.5915750915750916, 0.5915750915750916, 0.5970695970695971, 0.6007326007326007,
                 0.6007326007326007, 0.6043956043956044, 0.6043956043956044, 0.6117216117216118, 0.6117216117216118, 0.6135531135531136, 
                 0.6153846153846154, 0.6172161172161172, 0.6172161172161172, 0.6245421245421245, 0.6263736263736264, 0.6282051282051282,
                 0.63003663003663, 0.6318681318681318, 0.6373626373626373, 0.6391941391941391, 0.6410256410256411, 0.6428571428571429, 
                 0.6428571428571429, 0.6501831501831502, 0.6501831501831502, 0.6538461538461539, 0.6538461538461539, 0.6593406593406593,
                 0.663003663003663, 0.663003663003663, 0.6703296703296703, 0.6703296703296703, 0.6721611721611722, 0.673992673992674,
                 0.6794871794871795, 0.6813186813186813, 0.6831501831501832, 0.6868131868131868, 0.6868131868131868, 0.6923076923076923, 
                 0.6923076923076923, 0.6959706959706959, 0.6959706959706959, 0.7014652014652014, 0.7051282051282052, 0.7051282051282052, 
                 0.7124542124542125, 0.7124542124542125, 0.7161172161172161, 0.7161172161172161, 0.7216117216117216, 0.7234432234432234, 
                 0.7252747252747253, 0.7252747252747253, 0.7307692307692307, 0.7344322344322345, 0.7344322344322345, 0.7417582417582418, 
                 0.7417582417582418, 0.7454212454212454, 0.7472527472527473, 0.7472527472527473, 0.7545787545787546, 0.7545787545787546, 
                 0.7619047619047619, 0.7619047619047619, 0.7637362637362637, 0.7655677655677655, 0.7710622710622711, 0.7710622710622711, 
                 0.7728937728937729, 0.7802197802197802, 0.7802197802197802, 0.7838827838827839, 0.7838827838827839, 0.7912087912087912, 
                 0.793040293040293, 0.793040293040293, 0.7967032967032966, 0.7967032967032966, 0.8040293040293041, 0.8040293040293041, 
                 0.8058608058608059, 0.8076923076923077, 0.8131868131868132, 0.8168498168498168, 0.8168498168498168, 0.8223443223443223,
                 0.8223443223443223, 0.826007326007326, 0.826007326007326, 0.8296703296703297, 0.8315018315018315, 0.836996336996337, 
                 0.8388278388278388, 0.8388278388278388, 0.8424908424908425, 0.8424908424908425, 0.8498168498168498, 0.8498168498168498, 
                 0.8553113553113553, 0.8589743589743589, 0.8589743589743589, 0.8626373626373627, 0.8626373626373627, 0.8663003663003663, 
                 0.8663003663003663, 0.8717948717948718, 0.8736263736263736, 0.8754578754578755, 0.8809523809523809, 0.8827838827838828,
                 0.8846153846153846, 0.8846153846153846, 0.8882783882783882, 0.8882783882783882, 0.891941391941392, 0.8937728937728938, 
                 0.8956043956043956, 0.8974358974358975, 0.8974358974358975, 0.9010989010989011, 0.9010989010989011, 0.9084249084249084, 
                 0.9084249084249084, 0.9120879120879121, 0.9175824175824175, 0.9175824175824175, 0.9212454212454212, 0.9212454212454212, 
                 0.924908424908425, 0.924908424908425, 0.9267399267399268, 0.9267399267399268, 0.9304029304029304, 0.9340659340659341, 
                 0.9340659340659341, 0.9377289377289377, 0.9377289377289377, 0.9413919413919414, 0.9413919413919414, 0.9432234432234432, 
                 0.9468864468864469, 0.9468864468864469, 0.9468864468864469, 0.9468864468864469, 0.9505494505494505, 0.9505494505494505, 
                 0.9542124542124543, 0.9542124542124543, 0.9560439560439561, 0.9597069597069597, 0.9597069597069597, 0.9633699633699634, 
                 0.9633699633699634, 0.967032967032967, 0.967032967032967, 0.967032967032967, 0.967032967032967, 0.967032967032967, 
                 0.9706959706959707, 0.9706959706959707, 0.9725274725274725, 0.9725274725274725, 0.9725274725274725, 0.9743589743589743,
                 0.9761904761904762, 0.9761904761904762, 0.9761904761904762, 0.9798534798534798, 0.9798534798534798, 0.9798534798534798, 
                 0.9798534798534798, 0.9835164835164835, 0.9853479853479854, 0.9853479853479854, 0.9853479853479854, 0.9871794871794872, 
                 0.989010989010989, 0.989010989010989, 0.989010989010989, 0.989010989010989, 0.989010989010989, 0.9926739926739927,
                 0.9926739926739927, 0.9926739926739927, 0.9926739926739927, 0.9926739926739927, 0.9926739926739927, 0.9926739926739927, 
                 0.9945054945054945, 0.9963369963369964, 0.9963369963369964, 0.9963369963369964, 0.9963369963369964, 0.9963369963369964, 
                 0.9963369963369964, 0.9963369963369964, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 0.9963369963369964, 0.9963369963369964, 
                 1.0, 1.0, 1.0, 0.9963369963369964, 0.9963369963369964, 0.9963369963369964, 0.9963369963369964, 0.9963369963369964, 
                 0.9963369963369964, 0.9926739926739927, 0.9926739926739927, 0.9926739926739927, 0.989010989010989, 0.989010989010989, 
                 0.989010989010989, 0.989010989010989, 0.9871794871794872, 0.9853479853479854, 0.9853479853479854, 0.9835164835164835, 
                 0.9835164835164835, 0.9798534798534798, 0.9798534798534798, 0.9798534798534798, 0.9798534798534798, 0.9798534798534798,
                 0.978021978021978, 0.9761904761904762, 0.9761904761904762, 0.9761904761904762, 0.9725274725274725, 0.9725274725274725,
                 0.9706959706959707, 0.9688644688644689, 0.967032967032967, 0.9633699633699634, 0.9633699633699634, 0.9597069597069597, 
                 0.9597069597069597, 0.9597069597069597, 0.9597069597069597, 0.9560439560439561, 0.9560439560439561, 0.9542124542124543, 
                 0.9505494505494505, 0.9505494505494505, 0.9468864468864469, 0.9468864468864469, 0.9432234432234432, 0.9432234432234432, 
                 0.9377289377289377, 0.9340659340659341, 0.9340659340659341, 0.9304029304029304, 0.9304029304029304, 0.924908424908425, 
                 0.924908424908425, 0.9212454212454212, 0.9212454212454212, 0.9175824175824175, 0.9120879120879121, 0.9120879120879121, 
                 0.9084249084249084, 0.9084249084249084, 0.9010989010989011, 0.8974358974358975, 0.8974358974358975, 0.891941391941392, 
                 0.891941391941392, 0.8882783882783882, 0.8882783882783882, 0.8827838827838828, 0.8809523809523809, 0.8791208791208791, 
                 0.8791208791208791, 0.8681318681318682, 0.8626373626373627, 0.8626373626373627, 0.8589743589743589, 0.8589743589743589,
                 0.8516483516483516, 0.8516483516483516, 0.8461538461538461, 0.836996336996337, 0.836996336996337, 0.8333333333333334, 
                 0.8333333333333334, 0.826007326007326, 0.826007326007326, 0.8205128205128205, 0.8205128205128205, 0.8131868131868132, 
                 0.8076923076923077, 0.8058608058608059, 0.8003663003663004, 0.8003663003663004, 0.793040293040293, 0.793040293040293,
                 0.7875457875457875, 0.7765567765567766, 0.7765567765567766, 0.7673992673992674, 0.7673992673992674, 0.7619047619047619,
                 0.7619047619047619, 0.7545787545787546, 0.7545787545787546, 0.7472527472527473, 0.7417582417582418, 0.7417582417582418,
                 0.7326007326007326, 0.7326007326007326, 0.7252747252747253, 0.7252747252747253, 0.717948717948718, 0.717948717948718, 
                 0.7087912087912088, 0.7032967032967034, 0.7014652014652014, 0.6959706959706959, 0.6959706959706959, 0.6868131868131868,
                 0.6868131868131868, 0.6758241758241759, 0.6703296703296703, 0.6703296703296703, 0.6611721611721612, 0.6593406593406593, 
                 0.6501831501831502, 0.6501831501831502, 0.6428571428571429, 0.6428571428571429, 0.6373626373626373, 0.6282051282051282, 
                 0.6263736263736264, 0.6208791208791209, 0.6208791208791209, 0.6117216117216118, 0.6117216117216118, 0.6007326007326007, 
                 0.6007326007326007, 0.5952380952380952, 0.5860805860805861, 0.5842490842490843, 0.5787545787545788, 0.5787545787545788, 
                 0.5677655677655677, 0.5677655677655677, 0.5622710622710623, 0.5531135531135531, 0.5512820512820513, 0.5457875457875457, 
                 0.5457875457875457, 0.5366300366300366, 0.5366300366300366, 0.5256410256410257, 0.5256410256410257, 0.5201465201465202, 
                 0.510989010989011, 0.5091575091575091, 0.5036630036630036, 0.5036630036630036, 0.4926739926739927, 0.4926739926739927, 
                 0.4835164835164835, 0.47802197802197804, 0.47802197802197804, 0.4706959706959707, 0.4706959706959707, 0.46153846153846156,
                 0.46153846153846156, 0.45054945054945056, 0.44505494505494503, 0.44505494505494503, 0.43772893772893773, 
                 0.43772893772893773, 0.42857142857142855, 0.42857142857142855, 0.42124542124542125, 0.42124542124542125, 
                 0.41208791208791207, 0.4065934065934066, 0.40476190476190477, 0.3956043956043956, 0.3956043956043956, 
                 0.3901098901098901, 0.3882783882783883, 0.38278388278388276, 0.38278388278388276, 0.37545787545787546,
                 0.3663003663003663, 0.3663003663003663, 0.3608058608058608, 0.358974358974359, 0.3534798534798535, 
                 0.3534798534798535, 0.3424908424908425, 0.3424908424908425, 0.336996336996337, 0.3315018315018315, 
                 0.32967032967032966, 0.3241758241758242, 0.3241758241758242, 0.31684981684981683, 0.31135531135531136, 
                 0.30952380952380953, 0.3058608058608059, 0.304029304029304, 0.29853479853479853, 0.29853479853479853, 
                 0.29120879120879123, 0.29120879120879123, 0.2838827838827839, 0.2838827838827839, 0.2783882783882784, 
                 0.27289377289377287, 0.27106227106227104, 0.26556776556776557, 0.26556776556776557, 0.25824175824175827,
                 0.25824175824175827, 0.25274725274725274, 0.24725274725274726, 0.2454212454212454, 0.24175824175824176, 
                 0.24175824175824176, 0.23626373626373626, 0.23626373626373626, 0.22893772893772893, 0.22893772893772893,
                 0.22344322344322345, 0.21978021978021978, 0.21978021978021978, 0.21245421245421245, 0.21245421245421245, 
                 0.20695970695970695, 0.20695970695970695, 0.2032967032967033, 0.2032967032967033, 0.19597069597069597,
                 0.19413919413919414, 0.19230769230769232, 0.18681318681318682, 0.18681318681318682, 0.18315018315018314,
                 0.18315018315018314, 0.17765567765567766, 0.17399267399267399, 0.17399267399267399, 0.16666666666666666, 
                 0.16666666666666666, 0.16483516483516483, 0.163003663003663, 0.1575091575091575, 0.15384615384615385, 
                 0.15384615384615385, 0.152014652014652, 0.15018315018315018, 0.14835164835164835, 0.14835164835164835, 
                 0.14102564102564102, 0.14102564102564102, 0.13736263736263737, 0.13553113553113552, 0.1336996336996337,
                 0.13186813186813187, 0.13186813186813187, 0.1282051282051282, 0.1282051282051282, 0.12454212454212454, 
                 0.12454212454212454, 0.11904761904761904, 0.11538461538461539, 0.11538461538461539, 0.11172161172161173,
                 0.11172161172161173, 0.10805860805860806, 0.10805860805860806, 0.1043956043956044, 0.10256410256410256, 
                 0.10073260073260074, 0.0989010989010989, 0.0989010989010989, 0.09523809523809523, 0.09523809523809523, 
                 0.09157509157509157, 0.09157509157509157, 0.09157509157509157, 0.08974358974358974, 0.08791208791208792,
                 0.08608058608058608, 0.08608058608058608, 0.08241758241758242, 0.08241758241758242, 0.07875457875457875,
                 0.07875457875457875, 0.07509157509157509, 0.07326007326007326, 0.07326007326007326, 0.07326007326007326, 
                 0.07326007326007326, 0.0695970695970696, 0.0695970695970696, 0.06593406593406594, 0.0641025641025641, 
                 0.06227106227106227, 0.06227106227106227, 0.06227106227106227, 0.06043956043956044, 0.05860805860805861,
                 0.056776556776556776, 0.056776556776556776, 0.056776556776556776, 0.05311355311355311, 0.05311355311355311, 
                 0.05311355311355311, 0.05311355311355311, 0.04945054945054945, 0.04945054945054945, 0.04945054945054945,
                 0.047619047619047616, 0.045787545787545784, 0.04395604395604396, 0.04395604395604396, 0.04395604395604396, 
                 0.04395604395604396, 0.040293040293040296, 0.040293040293040296, 0.040293040293040296, 0.040293040293040296,
                 0.040293040293040296, 0.03663003663003663, 0.03663003663003663, 0.03663003663003663, 0.03663003663003663, 
                 0.03296703296703297, 0.03296703296703297, 0.03296703296703297, 0.031135531135531136, 0.029304029304029304, 
                 0.029304029304029304, 0.029304029304029304, 0.029304029304029304, 0.029304029304029304, 0.027472527472527472,
                 0.027472527472527472, 0.027472527472527472, 0.027472527472527472, 0.027472527472527472, 0.023809523809523808,
                 0.023809523809523808, 0.023809523809523808, 0.023809523809523808, 0.023809523809523808, 0.020146520146520148, 
                 0.020146520146520148, 0.020146520146520148, 0.020146520146520148, 0.020146520146520148, 0.020146520146520148, 
                 0.016483516483516484, 0.016483516483516484, 0.016483516483516484, 0.016483516483516484, 0.016483516483516484, 
                 0.016483516483516484, 0.016483516483516484, 0.016483516483516484, 0.016483516483516484, 0.01282051282051282, 
                 0.01282051282051282, 0.01282051282051282, 0.01282051282051282, 0.01282051282051282, 0.01282051282051282,
                 0.01282051282051282, 0.01282051282051282, 0.01098901098901099, 0.01098901098901099, 0.01098901098901099, 
                 0.01098901098901099, 0.01098901098901099, 0.01098901098901099, 0.01098901098901099, 0.01098901098901099, 
                 0.01098901098901099, 0.01098901098901099, 0.01098901098901099, 0.007326007326007326, 0.007326007326007326,
                 0.007326007326007326, 0.007326007326007326, 0.007326007326007326, 0.007326007326007326, 0.007326007326007326,
                 0.007326007326007326, 0.007326007326007326, 0.007326007326007326, 0.007326007326007326]

print("THORLABS SLD850S-A20W 提取的 波长vs强度 谱域特性曲线 (800-900nm)")
plt.figure()
index = [800+ (900-800)*i/(len(specdata_wlen)-1) for i in range(len(specdata_wlen))]
plt.plot(index,specdata_wlen)
freqindex = [i*376/511 for i in range(512)] #从0Hz ~ 400e+6MHz
ifreqindex = np.array(freqindex,dtype=int)
lightspeed = 3e+5 #光速 3e+5 米/毫秒
specdata_freq = [0.0]*1024
for i in range(1,512):
    t = freqindex[i]
    t = lightspeed / t #求波长
    #按波长查找幅值
    if t < 800:
        continue
    if t > 900:
        continue
    t = (t - 800) / (900 - 800)
    t = t * (len(specdata_wlen) - 1)
    j = int(t)
    d = t - j
    if j + 1 >= len(specdata_wlen):
        j = len(specdata_wlen) - 1
        d = 0
    t = specdata_wlen[j]*(1-d) + specdata_wlen[j + 1]*(d)
    specdata_freq[i] = t

print("转换为频谱图 x轴单位:1e+6 MHz")
plt.figure()
plt.plot(ifreqindex,specdata_freq[:len(specdata_freq)//2])

#让频谱图对称
for i in range(1,512):
    specdata_freq[1024-i] = specdata_freq[i]
#plt.figure()
#plt.plot(specdata_freq)
sldlightwave = np.fft.ifft(specdata_freq)
sldlightwave = np.concatenate((sldlightwave[-512:],sldlightwave[:512]))

sindex = [""]*1024
eversteptime = 1/(1024e+6) #每个ifft结果的时间间隔 单位us
for i in range(1024):
    t = (512 - i) * eversteptime
    t = t * lightspeed #单位 = us * m/ms = 0.001m = mm
    t = t * 10000 #转换为 0.1um
    sindex[i] = t
    #print(sindex[i])

print("理论上的光的波动形状图 x轴单位:0.1微米 0.1um")
#plt.figure()
plt.figure(figsize=(10,5))
plt.xlim(-150,150)
plt.plot(sindex,sldlightwave.real)

THORLABS SLD850S-A20W 提取的 波长vs强度 谱域特性曲线 (800-900nm)
转换为频谱图 x轴单位:1e+6 MHz
理论上的光的波动形状图 x轴单位:0.1微米 0.1um

在这里插入图片描述在这里插入图片描述在这里插入图片描述

如果TD-OCT光学系统在参考臂电机移动情况下光电二极管采样,应该会接收到下面模拟的波形(假设只有一个位置全反射的样品,反射的强度为100%)

print("相移后合成的波形")
srcdata = x
#srcdata = sldlightwave.real[256:768] #如果屏幕就使用正态函数乘以正弦波 的波形
phaseshift = 0
i = 0
start = -22
end = 22
cnt = end - start + 1 
plt.figure(figsize=(20,20))
amp = [0.0] * cnt
for phaseshift in range(start,end + 1,1):
    xshift = np.concatenate((srcdata[-phaseshift:], srcdata[:-phaseshift]))
    add = xshift + srcdata
    amp[i] = np.mean(np.sqrt(add**2))
    i = i + 1
    plt.subplot(cnt,1,i)
    plt.plot(add)
    plt.text(0,0,str(phaseshift))

#print("合成波形的强度图")
#plt.figure()
#plt.plot(amp)

相移后合成的波形
在这里插入图片描述

估计光电二极管转换为电信号应该是波形的强度(跟均根方值成比例吧),所以我估计可能光电二极管转换为电信号实际采样到波形会像下面模拟的这样(假设只有一个全反射的样品,反射的强度为20%)

#print("相移后合成的波形")
phaseshift = 0
i = 0
start = -128
end = 128
cnt = end - start + 1 
#plt.figure(figsize=(20,20))
amp = [0.0] * cnt
amp1 = [0.0] * cnt
for phaseshift in range(start,end + 1,1):
    xshift = np.concatenate((x[-phaseshift:], x[:-phaseshift]))
    add = 0.2 * xshift + srcdata
    amp[i] = np.mean(np.sqrt(add**2))
    amp1[i] = np.mean(np.sqrt(np.power(add[u+start//2 : u + end//2 +1],2)))
    i = i + 1
    #plt.subplot(cnt,1,i)
    #plt.plot(add)
    #plt.text(0,0,str(phaseshift))

print("合成波形的强度图(均根方值)")
print("左边相当于采样宽度很宽,右边相当于模拟实际光电二极管采样响应(一段时间内有效)")
plt.figure()
plt.subplot(1,2,1)
plt.plot(amp)
plt.subplot(1,2,2)
plt.plot(amp1)

合成波形的强度图(均根方值)
左边相当于采样宽度很宽,右边相当于模拟实际光电二极管采样响应(一段时间内有效)
在这里插入图片描述

从模拟结果来看,两个图在等光程(相移为0)时的波形振幅是最强的,这样的话跟上面找到的文章里提到的原理是吻合吧。这里我估计实际波形应该会像右图吧(光电二极管可能有转换周期之类,或者ADC有采样周期之类)。

问同事说用希尔伯特变换可以复原深度反射信息
先查资料找到一个个人认为比较有用的信息:
https://www.bilibili.com/read/cv2793412/

网页上,希尔伯特变换说是 1 t \frac{1}{t} t1傅里叶变换是 -j和j,我对 e − j w t t \frac{e^{-jwt}}{t} tejwt不会求积分,所以就不看公式了。希尔伯特变换说是对信号低频率部分相移 − π 2 -\frac{\pi}{2} 2π,高频部分相移 π 2 \frac{\pi}{2} 2π
希尔伯特变换用在信号解调方面的应用原理如下:
假设我们有一个低频信号a(t),用高频载波信号cos(wt)对其调制结果为 a ( t ) c o s ( w t ) a(t)cos(wt) a(t)cos(wt),如果我们能从 a ( t ) c o s ( w t ) a(t)cos(wt) a(t)cos(wt)经过滤波可以得到 a ( t ) s i n ( w t ) a(t)sin(wt) a(t)sin(wt)的话,就可以合成一个复数信号:
a ( t ) c o s ( w t ) + j a ( t ) s i n ( w t ) = a ( t ) ∗ ( c o s ( w t ) + j s i n ( w t ) ) = a ( t ) ∗ e j w t a(t)cos(wt)+ja(t)sin(wt) = a(t)*(cos(wt)+jsin(wt))=a(t)*e^{jwt} a(t)cos(wt)+ja(t)sin(wt)=a(t)(cos(wt)+jsin(wt))=a(t)ejwt
看公式就知到,要解调出原始信号a(t),只需要对复数信号 a ( t ) c o s ( w t ) + j a ( t ) s i n ( w t ) = a ( t ) ∗ e j w t a(t)cos(wt)+ja(t)sin(wt)=a(t)*e^{jwt} a(t)cos(wt)+ja(t)sin(wt)=a(t)ejwt 乘以 e − j w t e^{-jwt} ejwt 就行了。至于怎么样从调制信号 a ( t ) c o s ( w t ) a(t)cos(wt) a(t)cos(wt) 获得 a ( t ) s i n ( w t ) a(t)sin(wt) a(t)sin(wt) 应该就是希尔伯特变换显身手的地方了。
先来复习一下欧拉公式: e j x = c o s x + j s i n x e^{jx} = cosx + jsinx ejx=cosx+jsinx
再平复习一下共轭复数: ( x + j y ) ∗ ( x − j y ) = x 2 + y 2 (x+jy)*(x-jy) = x^{2}+y^{2} (x+jy)(xjy)=x2+y2
又复习一下离散傅里叶: X ( k ) = ∑ n = 0 N − 1 x n ∗ e − j 2 π k n N X(k) = \sum_{n=0}^{N-1}x_{n}*e^{\frac{-j2\pi kn}{N}} X(k)=n=0N1xneNj2πkn
离散傅里叶结果相对于 k = N 2 k=\frac{N}{2} k=2N共轭对称的,所以如果 k ′ = 1 , . . N 2 − 1 k' = 1,.. \frac{N}{2}-1 k=1,..2N1,有:
X ( k ′ ) = A m p 2 ( X ( N − k ′ ) ) X ( N − k ′ ) X(k') = \frac{Amp^2(X(N-k'))}{X(N-k')} X(k)=X(Nk)Amp2(X(Nk)) ,Amp()代表幅值

e − j 2 π k n N e^{\frac{-j2\pi kn}{N}} eNj2πkn W N k n W_N^{kn} WNkn
X ( k ) = ∑ n = 0 N − 1 x n ∗ e − j 2 π k n N = ∑ n = 0 N − 1 x n W N k n X(k) = \sum_{n=0}^{N-1}x_{n}*e^{\frac{-j2\pi kn}{N}} = \sum_{n=0}^{N-1}x_{n}W_N^{kn} X(k)=n=0N1xneNj2πkn=n=0N1xnWNkn
由于 W N − k n = W N − k n ∗ W N n N = W N ( N − k ) n W_N^{-kn}=W_N^{-kn} *W_N^{nN}=W_N^{(N-k)n} WNkn=WNknWNnN=WN(Nk)n可知,当 k > N 2 k>\frac{N}{2} k>2N时的X(k)是 k < N 2 k<\frac{N}{2} k<2N负频率共轭复数,所以:
X ( k ′ ) = A m p 2 ( X ( N − k ′ ) ) X ( N − k ′ ) = A m p 2 ( X ( k ′ ) ) X ( − k ′ ) X(k') = \frac{Amp^2(X(N-k'))}{X(N-k')}= \frac{Amp^2(X(k'))}{X(-k')} X(k)=X(Nk)Amp2(X(Nk))=X(k)Amp2(X(k))

将高频信号 c o s ( 2 π c n N ) cos(\frac{2\pi cn}{N}) cos(N2πcn)表示成复数的形式:
c o s ( 2 π c n N ) = e − j 2 π c n N + e j 2 π c n N 2 = W N c n + W N − c n 2 cos(\frac{2\pi cn}{N}) = \frac{e^{\frac{-j2\pi cn}{N}} + e^{\frac{j2\pi cn}{N}} }{2} = \frac{W_N^{cn}+W_N^{-cn}}{2} cos(N2πcn)=2eNj2πcn+eNj2πcn=2WNcn+WNcn

原始信号 x n x_n xn 经高频信号 c o s ( 2 π c n N ) cos(\frac{2\pi cn}{N}) cos(N2πcn) 调制后的信号为
x n ∗ W N c n + W N − c n 2 x_{n}*\frac{W_N^{cn}+W_N^{-cn}}{2} xn2WNcn+WNcn,它的傅里叶变换为:
X ′ ( k ) = ∑ n = 0 N − 1 x n ∗ W N c n + W N − c n 2 ∗ W N k n X'(k) = \sum_{n=0}^{N-1}x_{n}*\frac{W_N^{cn}+W_N^{-cn}}{2}*W_N^{kn} X(k)=n=0N1xn2WNcn+WNcnWNkn
对其进行化简
X ′ ( k ) = ∑ n = 0 N − 1 x n ∗ W N c n + W N − c n 2 ∗ W N k n = 1 2 ∑ n = 0 N − 1 x n ( W N c n + W N − c n ) W N k n = 1 2 ∑ n = 0 N − 1 x n ( W N ( k + c ) n + W N ( k − c ) n ) = 1 2 ∑ n = 0 N − 1 x n W N ( k + c ) n + 1 2 ∑ n = 0 N − 1 x n W N ( k − c ) n \begin{aligned} X'(k) &= \sum_{n=0}^{N-1}x_{n}*\frac{W_N^{cn}+W_N^{-cn}}{2}*W_N^{kn} \hspace{4cm} \\ &=\frac{1}{2} \sum_{n=0}^{N-1}x_{n}(W_N^{cn}+W_N^{-cn})W_N^{kn} \\ &=\frac{1}{2} \sum_{n=0}^{N-1}x_{n}(W_N^{(k+c)n}+W_N^{(k-c)n}) \\ &=\frac{1}{2} \sum_{n=0}^{N-1}x_{n}W_N^{(k+c)n} + \frac{1}{2} \sum_{n=0}^{N-1}x_{n}W_N^{(k-c)n} \end{aligned} X(k)=n=0N1xn2WNcn+WNcnWNkn=21n=0N1xn(WNcn+WNcn)WNkn=21n=0N1xn(WN(k+c)n+WN(kc)n)=21n=0N1xnWN(k+c)n+21n=0N1xnWN(kc)n

可以看到 被调制后的信号分成了 1 2 ∑ n = 0 N − 1 x n W N ( k + c ) n \frac{1}{2} \sum_{n=0}^{N-1}x_{n}W_N^{(k+c)n} 21n=0N1xnWN(k+c)n 1 2 ∑ n = 0 N − 1 x n W N ( k − c ) n \frac{1}{2} \sum_{n=0}^{N-1}x_{n}W_N^{(k-c)n} 21n=0N1xnWN(kc)n两部分

假设原始信号 x n x_n xn只是含有m个低频频率的信号(m小于 N 2 \frac{N}{2} 2N),那么它的傅里叶变换只含有 X ( − m + 1 ) , X ( − m + 2 ) . . . X ( − 1 ) , X ( 0 ) , X ( 1 ) , . . . X ( m − 1 ) X(-m+1),X(-m+2)...X(-1),X(0),X(1),...X(m-1) X(m+1),X(m+2)...X(1),X(0),X(1),...X(m1),其余频率X(m++或 - -m) 为0,

负频率 X ( − m + 1 ) , X ( − m + 2 ) . . . X ( − 1 ) X(-m+1),X(-m+2)...X(-1) X(m+1),X(m+2)...X(1)
相当于 X ( N − m + 1 ) , X ( N − m + 2 ) . . . X ( N − 1 ) X(N-m+1),X(N-m+2)...X(N-1) X(Nm+1),X(Nm+2)...X(N1)
换句话说 x n x_n xn只含有0,1,2…m 和 N-m+1,N-m+2,N-m+3…N-1两部分的频谱

假设高频信号 c o s ( 2 π c n N ) cos(\frac{2\pi cn}{N}) cos(N2πcn) 的c远大于m又小于 N 2 \frac{N}{2} 2N

对于 1 2 ∑ n = 0 N − 1 x n W N ( k + c ) n \frac{1}{2} \sum_{n=0}^{N-1}x_{n}W_N^{(k+c)n} 21n=0N1xnWN(k+c)n

可知 c ≤ k + c < ( N + N 2 ) c \le k+c<(N+\frac{N}{2}) ck+c<(N+2N)


c ≤ k + c ≤ N − m c \le k+c \le N-m ck+cNm 时( c > m c>m c>m),

1 2 ∑ n = 0 N − 1 x n W N ( k + c ) n = 1 2 X ( k + c ) = 0 \frac{1}{2} \sum_{n=0}^{N-1}x_{n}W_N^{(k+c)n} = \frac{1}{2}X(k+c)=0 21n=0N1xnWN(k+c)n=21X(k+c)=0


N − m + 1 ≤ k + c ≤ N − 1 N-m+1 \le k+c \le N-1 Nm+1k+cN1,

1 2 ∑ n = 0 N − 1 x n W N ( k + c ) n = 1 2 X ( k + c ) \frac{1}{2} \sum_{n=0}^{N-1}x_{n}W_N^{(k+c)n} = \frac{1}{2}X(k+c) 21n=0N1xnWN(k+c)n=21X(k+c)


N ≤ k + c ≤ N + m − 1 N \le k+c \le N+m-1 Nk+cN+m1 时( 0 ≤ k + c − N ≤ m − 1 0 \le k+c-N \le m-1 0k+cNm1),

1 2 ∑ n = 0 N − 1 x n W N ( k + c ) n = 1 2 ∑ n = 0 N − 1 x n W N ( k + c ) n W N − n N = 1 2 X ( k + c − N ) \frac{1}{2} \sum_{n=0}^{N-1}x_{n}W_N^{(k+c)n} =\frac{1}{2} \sum_{n=0}^{N-1}x_{n}W_N^{(k+c)n}W_N^{-nN} = \frac{1}{2}X(k+c-N) 21n=0N1xnWN(k+c)n=21n=0N1xnWN(k+c)nWNnN=21X(k+cN)


N + m ≤ k + c < N + N 2 N+m \le k+c < N+\frac{N}{2} N+mk+c<N+2N 时( m ≤ k + c − N < N 2 m \le k+c-N < \frac{N}{2} mk+cN<2N),

1 2 ∑ n = 0 N − 1 x n W N ( k + c ) n = 1 2 ∑ n = 0 N − 1 x n W N ( k + c ) n W N − n N = 1 2 X ( k + c − N ) = 0 \frac{1}{2} \sum_{n=0}^{N-1}x_{n}W_N^{(k+c)n} =\frac{1}{2} \sum_{n=0}^{N-1}x_{n}W_N^{(k+c)n}W_N^{-nN} = \frac{1}{2}X(k+c-N)=0 21n=0N1xnWN(k+c)n=21n=0N1xnWN(k+c)nWNnN=21X(k+cN)=0



总结来说,对于 1 2 ∑ n = 0 N − 1 x n W N ( k + c ) n \frac{1}{2} \sum_{n=0}^{N-1}x_{n}W_N^{(k+c)n} 21n=0N1xnWN(k+c)n

当( N − m + 1 ≤ k + c ≤ N − 1 N-m+1 \le k+c \le N-1 Nm+1k+cN1)及( N ≤ k + c ≤ N + m − 1 N \le k+c \le N+m-1 Nk+cN+m1)时求和的值才不为零;而且,


当( N − c − m + 1 ≤ k ≤ N − c − 1 N-c-m+1 \le k \le N-c-1 Ncm+1kNc1),由于 m < < c < N 2 m<<c<\frac{N}{2} m<<c<2N时,这时的k总是在 ( N 2 , N − 1 ] (\frac{N}{2},N-1] (2N,N1]之间
1 2 ∑ n = 0 N − 1 x n W N ( k + c ) n \frac{1}{2} \sum_{n=0}^{N-1}x_{n}W_N^{(k+c)n} 21n=0N1xnWN(k+c)n 的值为原始信号 x n x_n xn傅里叶变换的X(N-m+1),X(N-m+2),X(N-m+3)…X(N-1)的二分之一


当( N − c ≤ k ≤ N − c + m − 1 N-c \le k \le N-c+m-1 NckNc+m1),由于 m < < c < N 2 m<<c<\frac{N}{2} m<<c<2N时,这时的k总是在 ( N 2 , N − 1 ] (\frac{N}{2},N-1] (2N,N1]之间
1 2 ∑ n = 0 N − 1 x n W N ( k + c ) n \frac{1}{2} \sum_{n=0}^{N-1}x_{n}W_N^{(k+c)n} 21n=0N1xnWN(k+c)n 的值为原始信号 x n x_n xn傅里叶变换的X(0),X(1),X(2)…X(m-1)的二分之一


对于 1 2 ∑ n = 0 N − 1 x n W N ( k + c ) n \frac{1}{2} \sum_{n=0}^{N-1}x_{n}W_N^{(k+c)n} 21n=0N1xnWN(k+c)n,它代表着被调制后的信号的大于 N 2 \frac{N}{2} 2N倍频的频谱(或者说是负频率),而且它的值是由原始信号的频谱合成的,按这样的排序
1 2 X ( N − m − 1 ) , 1 2 X ( N − m + 2 ) , 1 2 X ( N − m + 3 ) … 1 2 X ( N − 1 ) , 1 2 X ( 0 ) , 1 2 X ( 1 ) , 1 2 X ( 2 ) … X ( m − 1 ) \frac{1}{2}X(N-m-1),\frac{1}{2}X(N-m+2),\frac{1}{2}X(N-m+3)…\frac{1}{2}X(N-1),\frac{1}{2}X(0),\frac{1}{2}X(1),\frac{1}{2}X(2)…X(m-1) 21X(Nm1),21X(Nm+2),21X(Nm+3)21X(N1),21X(0),21X(1),21X(2)X(m1)
在调制后信号频谱的N-c频率处为中心进行放置



对于 1 2 ∑ n = 0 N − 1 x n W N ( k − c ) n \frac{1}{2} \sum_{n=0}^{N-1}x_{n}W_N^{(k-c)n} 21n=0N1xnWN(kc)n

可知 − c ≤ k − c < ( N − N 2 ) -c \le k-c<(N-\frac{N}{2}) ckc<(N2N)


− c ≤ k − c ≤ − m -c \le k-c \le -m ckcm 时( N − c ≤ N + k − c ≤ N − m N-c \le N+k-c \le N-m NcN+kcNm)( c > m , c < N 2 c>m,c<\frac{N}{2} c>m,c<2N),

1 2 ∑ n = 0 N − 1 x n W N ( k − c ) n = 1 2 X ( k − c ) = 0 \frac{1}{2} \sum_{n=0}^{N-1}x_{n}W_N^{(k-c)n} = \frac{1}{2}X(k-c)=0 21n=0N1xnWN(kc)n=21X(kc)=0


− m + 1 ≤ k − c < 0 -m+1 \le k-c < 0 m+1kc<0 时( N − m + 1 ≤ N + k − c < N N-m+1 \le N+k-c < N Nm+1N+kc<N)( c > m , c < N 2 c>m,c<\frac{N}{2} c>m,c<2N),

1 2 ∑ n = 0 N − 1 x n W N ( k − c ) n = 1 2 X ( k − c ) \frac{1}{2} \sum_{n=0}^{N-1}x_{n}W_N^{(k-c)n} = \frac{1}{2}X(k-c) 21n=0N1xnWN(kc)n=21X(kc)

对应序列 1 2 X ( N − m + 1 ) , 1 2 X ( N − m + 2 ) . . . 1 2 X ( N − 1 ) \frac{1}{2}X(N-m+1),\frac{1}{2}X(N-m+2)...\frac{1}{2}X(N-1) 21X(Nm+1),21X(Nm+2)...21X(N1)


0 ≤ k − c < m 0 \le k-c < m 0kc<m 时( c > m , c < N 2 c>m,c<\frac{N}{2} c>m,c<2N),

1 2 ∑ n = 0 N − 1 x n W N ( k − c ) n = 1 2 X ( k − c ) \frac{1}{2} \sum_{n=0}^{N-1}x_{n}W_N^{(k-c)n} = \frac{1}{2}X(k-c) 21n=0N1xnWN(kc)n=21X(kc)

对应序列 1 2 X ( 0 ) , 1 2 X ( 1 ) . . . 1 2 X ( m − 1 ) \frac{1}{2}X(0),\frac{1}{2}X(1)...\frac{1}{2}X(m-1) 21X(0),21X(1)...21X(m1)


m ≤ k − c < N − N 2 m \le k-c < N-\frac{N}{2} mkc<N2N 时( c > m , c < N 2 c>m,c<\frac{N}{2} c>m,c<2N),

1 2 ∑ n = 0 N − 1 x n W N ( k − c ) n = 1 2 X ( k − c ) = 0 \frac{1}{2} \sum_{n=0}^{N-1}x_{n}W_N^{(k-c)n} = \frac{1}{2}X(k-c)=0 21n=0N1xnWN(kc)n=21X(kc)=0



总结来说,对于 1 2 ∑ n = 0 N − 1 x n W N ( k − c ) n \frac{1}{2} \sum_{n=0}^{N-1}x_{n}W_N^{(k-c)n} 21n=0N1xnWN(kc)n

当( − m + 1 ≤ k − c < 0 -m+1 \le k-c < 0 m+1kc<0)及( 0 ≤ k − c < m 0 \le k-c < m 0kc<m)时求和的值才不为零;
也就是说( c − m + 1 ≤ k < c c-m+1 \le k < c cm+1k<c)及( c ≤ k < c + m c \le k < c+m ck<c+m)时,才有与原始信号频谱相关的值。

对于 1 2 ∑ n = 0 N − 1 x n W N ( k − c ) n \frac{1}{2} \sum_{n=0}^{N-1}x_{n}W_N^{(k-c)n} 21n=0N1xnWN(kc)n,它代表着被调制后的信号的小于 N 2 \frac{N}{2} 2N倍频的频谱,而且它的值是由原始信号的频谱合成的,按这样的排序
1 2 X ( N − m + 1 ) , 1 2 X ( N − m + 2 ) , 1 2 X ( N − m + 3 ) … 1 2 X ( N − 1 ) , 1 2 X ( 0 ) , 1 2 X ( 1 ) , 1 2 X ( 2 ) … X ( m − 1 ) \frac{1}{2}X(N-m+1),\frac{1}{2}X(N-m+2),\frac{1}{2}X(N-m+3)…\frac{1}{2}X(N-1),\frac{1}{2}X(0),\frac{1}{2}X(1),\frac{1}{2}X(2)…X(m-1) 21X(Nm+1),21X(Nm+2),21X(Nm+3)21X(N1),21X(0),21X(1),21X(2)X(m1)
在调制后信号频谱的c频率处为中心进行放置




所以对于调制后的信号

其傅里叶变换可化为
X ′ ( k ) = ∑ n = 0 N − 1 x n ∗ W N c n + W N − c n 2 ∗ W N k n = 1 2 ∑ n = 0 N − 1 x n W N ( k + c ) n + 1 2 ∑ n = 0 N − 1 x n W N ( k − c ) n \begin{aligned} X'(k) &= \sum_{n=0}^{N-1}x_{n}*\frac{W_N^{cn}+W_N^{-cn}}{2}*W_N^{kn} \hspace{4cm} \\ &=\frac{1}{2} \sum_{n=0}^{N-1}x_{n}W_N^{(k+c)n} + \frac{1}{2} \sum_{n=0}^{N-1}x_{n}W_N^{(k-c)n} \end{aligned} X(k)=n=0N1xn2WNcn+WNcnWNkn=21n=0N1xnWN(k+c)n+21n=0N1xnWN(kc)n
1 2 ∑ n = 0 N − 1 x n W N ( k − c ) n \frac{1}{2} \sum_{n=0}^{N-1}x_{n}W_N^{(k-c)n} 21n=0N1xnWN(kc)n代表着正半频率频谱, 1 2 ∑ n = 0 N − 1 x n W N ( k + c ) n \frac{1}{2} \sum_{n=0}^{N-1}x_{n}W_N^{(k+c)n} 21n=0N1xnWN(k+c)n代表着负半频率频谱(或者说是大于 N 2 \frac{N}{2} 2N倍频的频谱)
希尔伯特变换相当于正半频谱乘以-j,负半频谱乘以j,见下面推导:
H ( X ′ ( k ) ) = j 1 2 ∑ n = 0 N − 1 x n W N ( k + c ) n − j 1 2 ∑ n = 0 N − 1 x n W N ( k − c ) n = j 1 2 ∑ n = 0 N − 1 x n W N k n W N c n − j 1 2 ∑ n = 0 N − 1 x n W N k n W N − c n = ∑ n = 0 N − 1 x n j W N c n − j W N − c n 2 W N k n \begin{aligned} H(X'(k)) &=j\frac{1}{2} \sum_{n=0}^{N-1}x_{n}W_N^{(k+c)n} -j \frac{1}{2} \sum_{n=0}^{N-1}x_{n}W_N^{(k-c)n} \hspace{4cm} \\ &=j\frac{1}{2} \sum_{n=0}^{N-1}x_{n}W_N^{kn}W_N^{cn} -j \frac{1}{2} \sum_{n=0}^{N-1}x_{n}W_N^{kn}W_N^{-cn} \\ &= \sum_{n=0}^{N-1}x_{n}\frac{jW_N^{cn}-jW_N^{-cn}}{2}W_N^{kn} \end{aligned} H(X(k))=j21n=0N1xnWN(k+c)nj21n=0N1xnWN(kc)n=j21n=0N1xnWNknWNcnj21n=0N1xnWNknWNcn=n=0N1xn2jWNcnjWNcnWNkn
由上面的推度可知道
调制后的信号,傅里叶变换->希尔伯特变换->傅里叶逆变换 得到的是
x n j W N c n − j W N − c n 2 x_{n}\frac{jW_N^{cn}-jW_N^{-cn}}{2} xn2jWNcnjWNcn

j W N c n − j W N − c n 2 = j e − j 2 π c n N − j e j 2 π c n N 2 = 1 2 ( j ( c o s ( 2 π c n N ) − j s i n ( 2 π c n N ) − c o s ( 2 π c n N ) − j s i n ( 2 π c n N ) ) ) = 1 2 ( j ( − 2 ∗ j ∗ s i n ( 2 π c n N ) ) ) = s i n ( 2 π c n N ) \begin{aligned} \frac{jW_N^{cn}-jW_N^{-cn}}{2}&=\frac{je^{\frac{-j2\pi cn}{N}}-je^{\frac{j2\pi cn}{N}}}{2} \hspace{4cm}\\ &=\frac{1}{2}(j(cos(\frac{2\pi cn}{N})-jsin(\frac{2\pi cn}{N}) -cos(\frac{2\pi cn}{N})-jsin(\frac{2\pi cn}{N}) )) \\ &=\frac{1}{2}(j(-2*j*sin(\frac{2\pi cn}{N}))) \\ &=sin(\frac{2\pi cn}{N}) \end{aligned} 2jWNcnjWNcn=2jeNj2πcnjeNj2πcn=21(j(cos(N2πcn)jsin(N2πcn)cos(N2πcn)jsin(N2πcn)))=21(j(2jsin(N2πcn)))=sin(N2πcn)
也就是
x n j W N c n − j W N − c n 2 = x n s i n ( 2 π c n N ) x_{n}\frac{jW_N^{cn}-jW_N^{-cn}}{2}=x_{n}sin(\frac{2\pi cn}{N}) xn2jWNcnjWNcn=xnsin(N2πcn)




x n x_{n} xn通过高频载波信号 c o s ( 2 π c n N ) cos(\frac{2\pi cn}{N}) cos(N2πcn)调制后的信号 x n c o s ( 2 π c n N ) x_{n}cos(\frac{2\pi cn}{N}) xncos(N2πcn)进行解调方法:

  • x n c o s ( 2 π c n N ) x_{n}cos(\frac{2\pi cn}{N}) xncos(N2πcn)进行傅里叶变换得到 X ( k ) X(k) X(k)
  • k < N 2 k<\frac{N}{2} k<2N, X ′ ( k ) = − j X ( k ) X'(k)=-jX(k) X(k)=jX(k)
  • k > N 2 k>\frac{N}{2} k>2N, X ′ ( k ) = j X ( k ) X'(k)=jX(k) X(k)=jX(k)
  • X ′ ( k ) X'(k) X(k)进行傅里叶逆变换得到解析信号 x n s i n ( 2 π c n N ) x_{n}sin(\frac{2\pi cn}{N}) xnsin(N2πcn)
  • 将调制后的信号和解析信号合成复数信号 x n ′ = x n c o s ( 2 π c n N ) + j x n s i n ( 2 π c n N ) x'_{n} = x_{n}cos(\frac{2\pi cn}{N})+jx_{n}sin(\frac{2\pi cn}{N}) xn=xncos(N2πcn)+jxnsin(N2πcn)
  • 对合成的复数信号乘以 e − j 2 π c n N e^{\frac{-j2\pi cn}{N}} eNj2πcn,相乘后的结果实部原始信号 x n x_n xn
  • 如果不知道高频调制信号的频率c,可以对复数信号每个点求复数的幅值就是原始信号的幅值了

下面用python进行这个过程试验

a = [0] * 512
c = [0] * 512
data = [0] * 512
for i in range(512):
    #a[i] = i / 511
    sigma = 10
    u = 256
    a[i] = 1/(math.sqrt(2 * math.pi)*sigma) * math.exp(-(i-u)**2 / (2 * sigma**2))
    #a[i] = 1.0 - math.cos(2 * math.pi * i /511)
    c[i] = math.cos(2 * math.pi * 100 * i /511)
    data[i] = a[i] * c[i]

datafft = np.fft.fft(data)
    
plt.figure(figsize=(20,10))
plt.subplot(311)
plt.plot(np.abs(np.fft.fft(a)))

plt.subplot(312)
plt.plot(np.abs(np.fft.fft(c)))

plt.subplot(313)
plt.plot(np.abs(np.fft.fft(data)))


for i in range(1,256): # 希尔伯特变换 Hilbert transform
    datafft[i] = datafft[i] * (0-1j)
    datafft[i+256] = datafft[i+256] * (0+1j)
datashift = np.fft.ifft(datafft) # 希尔伯特变换 Hilbert transform


#print(datashift)
plt.figure(figsize=(20,5))
plt.plot(data)
plt.plot(datashift.real)

datacomplex = [0+0j]*512
for i in range(512):  #生成解释信号波型
    datacomplex[i] = data[i] + (0+1j) * datashift[i]

dataz = [0+0j]*512
for i in range(0,512):  #乘以e^-jwt 获取原始被调制信号
    dataz[i] = datacomplex[i] * cmath.exp((0-1j) * 2 * math.pi * 100 * i/511) 

plt.figure(figsize=(20,10))
plt.subplot(211)
plt.plot(a)
plt.subplot(212)
plt.plot(a)
plt.plot(data)
plt.plot(np.real(dataz)) #这个是在知道调制信号的频率的情况下才有用的
#plt.plot(np.abs(datacomplex)) #不知道调制信号,且只要求得到原始信号的幅值

在这里插入图片描述在这里插入图片描述在这里插入图片描述
可以看到代码执行后的结果与公式推导是吻合的。

用希尔伯特变换的方法对TD-OCT模拟的信号解调

#orgsignal = data
orgsignal = amp
#print(np.mean(orgsignal))
submean = np.array(orgsignal,dtype=np.float) - np.mean(orgsignal)
plt.plot(submean)

#进行 hilbert变换
hilbertfft = np.fft.fft(submean) 
hcnt = len(hilbertfft)//2

for i in range(hcnt):
    hilbertfft[i] = 0-1j * hilbertfft[i]
    hilbertfft[i+hcnt] = 0+1j * hilbertfft[i+hcnt]

#hilbertfft[0] = 0 #不要零频
#hilbertfft[hcnt] = 0

hilbertsignal = np.fft.ifft(hilbertfft)
plt.plot(hilbertsignal.real)

complexsignal = submean +  (0+1j * hilbertsignal)
decodesignal = np.abs(complexsignal)

plt.figure()
plt.plot(decodesignal)

#如像解调出来的波形并不理想
#进行均值滤波 效果好一点
cnvcnt = 7
cnvkernel = [1/cnvcnt]*cnvcnt
cnvsignal = np.convolve(decodesignal,cnvkernel,'same') #valid same full

plt.figure()
plt.plot(cnvsignal)

#其实好像直接用原信号(没有负数的信号)进行平均滤波效果更好  但对于实际TD-OCT采集的信号是不是也不好说
cnvsignal1 = np.convolve(orgsignal,cnvkernel,'same') #valid same full
plt.figure()
plt.plot(0.5-cnvsignal1)

在这里插入图片描述在这里插入图片描述
在这里插入图片描述在这里插入图片描述

希尔伯特变换 再深入

对于调制后的信号 x n ′ = x n ∗ c o s ( 2 π c n N ) x'_n = x_n * cos(\frac{2\pi cn}{N}) xn=xncos(N2πcn)

假设原始信号 x n x_n xn只是含有m个低频频率的信号(m小于 N 2 \frac{N}{2} 2N),那么它的傅里叶变换只含有 X ( − m + 1 ) , X ( − m + 2 ) . . . X ( − 1 ) , X ( 0 ) , X ( 1 ) , . . . X ( m − 1 ) X(-m+1),X(-m+2)...X(-1),X(0),X(1),...X(m-1) X(m+1),X(m+2)...X(1),X(0),X(1),...X(m1),其余频率X(m++或 - -m) 为0,

负频率 X ( − m + 1 ) , X ( − m + 2 ) . . . X ( − 1 ) X(-m+1),X(-m+2)...X(-1) X(m+1),X(m+2)...X(1)
相当于 X ( N − m + 1 ) , X ( N − m + 2 ) . . . X ( N − 1 ) X(N-m+1),X(N-m+2)...X(N-1) X(Nm+1),X(Nm+2)...X(N1)
换句话说 x n x_n xn只含有0,1,2…m 和 N-m+1,N-m+2,N-m+3…N-1两部分的频谱
假设高频载波信号 c o s ( 2 π c n N ) cos(\frac{2\pi cn}{N}) cos(N2πcn) 的c大于m又小于 N 2 \frac{N}{2} 2N

调制后的信号傅里叶变换可化为
X ′ ( k ) = ∑ n = 0 N − 1 x n ∗ W N c n + W N − c n 2 ∗ W N k n = 1 2 ∑ n = 0 N − 1 x n W N ( k + c ) n + 1 2 ∑ n = 0 N − 1 x n W N ( k − c ) n \begin{aligned} X'(k) &= \sum_{n=0}^{N-1}x_{n}*\frac{W_N^{cn}+W_N^{-cn}}{2}*W_N^{kn} \hspace{4cm} \\ &=\frac{1}{2} \sum_{n=0}^{N-1}x_{n}W_N^{(k+c)n} + \frac{1}{2} \sum_{n=0}^{N-1}x_{n}W_N^{(k-c)n} \end{aligned} X(k)=n=0N1xn2WNcn+WNcnWNkn=21n=0N1xnWN(k+c)n+21n=0N1xnWN(kc)n
1 2 ∑ n = 0 N − 1 x n W N ( k − c ) n \frac{1}{2} \sum_{n=0}^{N-1}x_{n}W_N^{(k-c)n} 21n=0N1xnWN(kc)n代表着正半频率频谱, 1 2 ∑ n = 0 N − 1 x n W N ( k + c ) n \frac{1}{2} \sum_{n=0}^{N-1}x_{n}W_N^{(k+c)n} 21n=0N1xnWN(k+c)n代表着负半频率频谱(或者说是大于 N 2 \frac{N}{2} 2N倍频的频谱)
希尔伯特变换相当于正半频谱乘以-j,负半频谱乘以j,
H i l b e r t ( x n ′ ) = H i l b e r t ( x n ∗ c o s ( 2 π c n N ) ) = x n c o s ( 2 π c n N − π 2 ) = x n s i n ( 2 π c n N ) Hilbert(x'_n) = Hilbert(x_n * cos(\frac{2\pi cn}{N}))=x_n cos(\frac{2\pi cn}{N}-\frac{\pi}{2})=x_n sin(\frac{2\pi cn}{N}) Hilbert(xn)=Hilbert(xncos(N2πcn))=xncos(N2πcn2π)=xnsin(N2πcn)

如果载波信号变为 C ∗ c o s ( 2 π c n N + ψ ) C * cos(\frac{2\pi cn}{N} + \psi) Ccos(N2πcn+ψ) 其中C为幅值, c相当于频率, ψ \psi ψ为相位
代入上面的公式依然有:
X ′ ( k ) = ∑ n = 0 N − 1 x n ∗ C ∗ W N ( c n + N ψ / 2 π ) + W N − ( c n + N ψ / 2 π ) 2 ∗ W N k n = 1 2 ∑ n = 0 N − 1 x n ∗ C ∗ ( W N c n e − j ψ + W N − c n e j ψ ) ∗ W N k n = 1 2 ∑ n = 0 N − 1 x n C W N ( k + c ) n e − j ψ + 1 2 ∑ n = 0 N − 1 x n C W N ( k − c ) n e j ψ = 1 2 ( c o s ψ − j s i n ψ ) ∑ n = 0 N − 1 x n C W N ( k + c ) n + 1 2 ( c o s ψ + j s i n ψ ) ∑ n = 0 N − 1 x n C W N ( k − c ) n = 1 2 ( ∑ n = 0 N − 1 x n c o s ψ C W N ( k + c ) n − j ∑ n = 0 N − 1 x n s i n ψ C W N ( k + c ) n ) + 1 2 ( ∑ n = 0 N − 1 c o s ψ x n C W N ( k − c ) n + j ∑ n = 0 N − 1 x n s i n ψ C W N ( k − c ) n ) \begin{aligned} X'(k) &= \sum_{n=0}^{N-1}x_{n}*C*\frac{W_N^{(cn+N\psi/2\pi)}+W_N^{-(cn+N\psi/2\pi)}}{2}*W_N^{kn} \hspace{4cm} \\ &=\frac{1}{2} \sum_{n=0}^{N-1}x_{n}*C*(W_N^{cn}e^{-j\psi}+W_N^{-cn}e^{j\psi})*W_N^{kn} \\ &=\frac{1}{2} \sum_{n=0}^{N-1}x_{n}CW_N^{(k+c)n}e^{-j\psi} + \frac{1}{2} \sum_{n=0}^{N-1}x_{n}CW_N^{(k-c)n}e^{j\psi} \\ &=\frac{1}{2} (cos\psi-jsin\psi)\sum_{n=0}^{N-1}x_{n}CW_N^{(k+c)n} + \frac{1}{2}(cos\psi+jsin\psi) \sum_{n=0}^{N-1}x_{n}CW_N^{(k-c)n}\\ &=\frac{1}{2}\left(\sum_{n=0}^{N-1} x_{n} cos\psi CW_N^{(k+c)n} - j \sum_{n=0}^{N-1}x_{n} sin\psi CW_N^{(k+c)n} \right) +\frac{1}{2}\left(\sum_{n=0}^{N-1} cos \psi x_{n}CW_N^{(k-c)n} + j\sum_{n=0}^{N-1}x_{n} sin\psi CW_N^{(k-c)n} \right) \\ \end{aligned} X(k)=n=0N1xnC2WN(cn+Nψ/2π)+WN(cn+Nψ/2π)WNkn=21n=0N1xnC(WNcnejψ+WNcnejψ)WNkn=21n=0N1xnCWN(k+c)nejψ+21n=0N1xnCWN(kc)nejψ=21(cosψjsinψ)n=0N1xnCWN(k+c)n+21(cosψ+jsinψ)n=0N1xnCWN(kc)n=21(n=0N1xncosψCWN(k+c)njn=0N1xnsinψCWN(k+c)n)+21(n=0N1cosψxnCWN(kc)n+jn=0N1xnsinψCWN(kc)n)
x n x_n xn乘以常数,频谱是不会变化的,k+c 和 k-c的特性还是没有改变,所以 1 2 ∑ n = 0 N − 1 x n C W N ( k − c ) n e j ψ \frac{1}{2} \sum_{n=0}^{N-1}x_{n}CW_N^{(k-c)n}e^{j\psi} 21n=0N1xnCWN(kc)nejψ代表着正半频率频谱, 1 2 ∑ n = 0 N − 1 x n C W N ( k + c ) n e − j ψ \frac{1}{2} \sum_{n=0}^{N-1}x_{n}CW_N^{(k+c)n}e^{-j\psi} 21n=0N1xnCWN(k+c)nejψ代表着负半频率频谱(或者说是大于 N 2 \frac{N}{2} 2N倍频的频谱)

调制后的信号->希尔伯特->傅里叶变换 有:
D F T ( H i l b e r t ( x n ′ ) ) = j 1 2 ∑ n = 0 N − 1 x n C W N ( k + c ) n e − j ψ − j 1 2 ∑ n = 0 N − 1 x n C W N ( k − c ) n e j ψ = 1 2 ∑ n = 0 N − 1 x n ∗ C ∗ j ( W N c n + N ψ / 2 π − W N − ( c n + N ψ / 2 π ) ) ∗ W N k n = 1 2 ∑ n = 0 N − 1 x n ∗ C ∗ j ( − 2 j s i n ( 2 π c n + N ψ N ) ) ∗ W N k n = ∑ n = 0 N − 1 x n ∗ C ∗ s i n ( 2 π c n N + ψ ) ∗ W N k n \begin{aligned} DFT(Hilbert(x'_n)) &=j\frac{1}{2} \sum_{n=0}^{N-1}x_{n}CW_N^{(k+c)n}e^{-j\psi} -j \frac{1}{2} \sum_{n=0}^{N-1}x_{n}CW_N^{(k-c)n}e^{j\psi} \hspace{4cm} \\ &=\frac{1}{2} \sum_{n=0}^{N-1}x_{n}*C*j(W_N^{cn+N\psi/2\pi}-W_N^{-(cn+N\psi/2\pi)})*W_N^{kn} \\ &=\frac{1}{2} \sum_{n=0}^{N-1}x_{n}*C*j(-2jsin(\frac{2\pi cn + N\psi}{N}))*W_N^{kn} \\ &=\sum_{n=0}^{N-1}x_{n}*C*sin(\frac{2\pi cn }{N}+ \psi)*W_N^{kn} \end{aligned} DFT(Hilbert(xn))=j21n=0N1xnCWN(k+c)nejψj21n=0N1xnCWN(kc)nejψ=21n=0N1xnCj(WNcn+Nψ/2πWN(cn+Nψ/2π))WNkn=21n=0N1xnCj(2jsin(N2πcn+Nψ))WNkn=n=0N1xnCsin(N2πcn+ψ)WNkn
也就是说
H i l b e r t ( x n ′ ) = H i l b e r t ( x n ∗ C ∗ c o s ( 2 π c n N + ψ ) ) = x n ∗ C ∗ s i n ( 2 π c n N + ψ ) = x n ∗ C ∗ c o s ( 2 π c n N + ψ − π 2 ) \begin{aligned} Hilbert(x'_n) &=Hilbert(x_{n}*C*cos(\frac{2\pi cn }{N}+ \psi)) \hspace{4cm} \\ &=x_{n}*C*sin(\frac{2\pi cn }{N}+ \psi) \\ &=x_{n}*C*cos(\frac{2\pi cn }{N}+ \psi - \frac{\pi}{2}) \\ \end{aligned} Hilbert(xn)=Hilbert(xnCcos(N2πcn+ψ))=xnCsin(N2πcn+ψ)=xnCcos(N2πcn+ψ2π)

如果载波信号由更多的余弦信号合成

∑ i = 0 I C i ∗ c o s ( 2 π c i n N + ψ i ) \sum_{i=0}^I C_i * cos(\frac{2\pi c_i n}{N} + \psi_i) i=0ICicos(N2πcin+ψi) 其中 C i Ci Ci为幅值, c i c_i ci相当于频率, ψ i \psi_i ψi为相位, c i c_i ci总是大于 x n x_n xn含有最大频率m,那么 x n x_n xn被调制后的信号 x n ′ x'_n xn可用下面式子表示:
x n ′ = x n ∗ ∑ i = 0 I C i ∗ c o s ( 2 π c i n N + ψ i ) \begin{aligned} x'_n &=x_{n}*\sum_{i=0}^I C_i * cos(\frac{2\pi c_i n}{N} + \psi_i) \hspace{4cm} \\ \end{aligned} xn=xni=0ICicos(N2πcin+ψi)
x n ′ x'_n xn作傅里叶变换可用下面式子表示:
X ′ ( k ) = 1 2 ∑ n = 0 N − 1 x n ∗ ( ∑ i = 0 i = I C i ∗ ( W N c i n + N ψ i / 2 π + W N − ( c i n + N ψ i / 2 π ) ) ) ∗ W N k n = 1 2 ∑ n = 0 N − 1 x n ( ∑ i = 0 i = I C i W N ( k + c i ) n e − j ψ i ) + 1 2 ∑ n = 0 N − 1 x n ( ∑ i = 0 i = I C i W N ( k − c i ) n e j ψ i ) = 1 2 ∑ n = 0 N − 1 ∑ i = 0 i = I x n C i W N ( k + c i ) n e − j ψ i + 1 2 ∑ n = 0 N − 1 ∑ i = 0 i = I x n C i W N ( k − c i ) n e j ψ i = 1 2 ∑ i = 0 i = I ∑ n = 0 N − 1 x n C i W N ( k + c i ) n e − j ψ i + 1 2 ∑ i = 0 i = I ∑ n = 0 N − 1 x n C i W N ( k − c i ) n e j ψ i \begin{aligned} X'(k) &=\frac{1}{2} \sum_{n=0}^{N-1}x_{n}* \left( \sum_{i=0}^{i=I} C_i*(W_N^{c_in+N\psi _i /2\pi}+W_N^{-(c_i n+N\psi _i /2\pi )}) \right)* W_N^{kn} \hspace{4cm} \\ &=\frac{1}{2} \sum_{n=0}^{N-1} x_{n} \left( \sum_{i=0}^{i=I} C_iW_N^{(k+c_i)n}e^{-j\psi_i} \right) + \frac{1}{2} \sum_{n=0}^{N-1} x_{n} \left( \sum_{i=0}^{i=I} C_iW_N^{(k-c_i)n}e^{j\psi_i} \right) \\ &=\frac{1}{2} \sum_{n=0}^{N-1}\sum_{i=0}^{i=I} x_{n} C_iW_N^{(k+c_i)n}e^{-j\psi_i} + \frac{1}{2} \sum_{n=0}^{N-1} \sum_{i=0}^{i=I}x_{n} C_iW_N^{(k-c_i)n}e^{j\psi_i} \\ &=\frac{1}{2} \sum_{i=0}^{i=I}\sum_{n=0}^{N-1} x_{n} C_iW_N^{(k+c_i)n}e^{-j\psi_i} + \frac{1}{2} \sum_{i=0}^{i=I}\sum_{n=0}^{N-1}x_{n} C_iW_N^{(k-c_i)n}e^{j\psi_i} \\ \end{aligned} X(k)=21n=0N1xn(i=0i=ICi(WNcin+Nψi/2π+WN(cin+Nψi/2π)))WNkn=21n=0N1xn(i=0i=ICiWN(k+ci)nejψi)+21n=0N1xn(i=0i=ICiWN(kci)nejψi)=21n=0N1i=0i=IxnCiWN(k+ci)nejψi+21n=0N1i=0i=IxnCiWN(kci)nejψi=21i=0i=In=0N1xnCiWN(k+ci)nejψi+21i=0i=In=0N1xnCiWN(kci)nejψi
这里的 k + c i k+c_i k+ci k − c i k-c_i kci与上面推导的特性也一样,分别对应着调制后信号的负频谱和正频谱。不同的是,原始信号被多个余弦信号合成后的信号调制,其傅里叶变换等于原始信号分别被多个余弦信号调制后傅里叶相加的结果 x n ′ x'_n xn调制后的信号希尔伯特变换再作傅里叶变换相当于:

D F T ( H i l b e r t ( x n ′ ) ) = j 1 2 ∑ i = 0 i = I ∑ n = 0 N − 1 x n C i W N ( k + c i ) n e − j ψ i − j 1 2 ∑ i = 0 i = I ∑ n = 0 N − 1 x n C i W N ( k − c i ) n e j ψ i = 1 2 ∑ n = 0 N − 1 x n ( ∑ i = 0 i = I j C i W N ( k + c i ) n e − j ψ i ) + 1 2 ∑ n = 0 N − 1 x n ( ∑ i = 0 i = I − j C i W N ( k − c i ) n e j ψ i ) = 1 2 ∑ n = 0 N − 1 x n ∗ ( ∑ i = 0 i = I C i ∗ j ( W N c i n + N ψ i / 2 π − W N − ( c i n + N ψ i / 2 π ) ) ) ∗ W N k n = 1 2 ∑ n = 0 N − 1 x n ∗ ( ∑ i = 0 i = I C i ∗ j ( − 2 j s i n ( 2 π c i n N + ψ i ) ) ) ∗ W N k n = ∑ n = 0 N − 1 x n ∗ ( ∑ i = 0 i = I C i ∗ s i n ( 2 π c i n N + ψ i ) ) ∗ W N k n \begin{aligned} DFT(Hilbert(x'_n)) &= j\frac{1}{2} \sum_{i=0}^{i=I}\sum_{n=0}^{N-1} x_{n} C_iW_N^{(k+c_i)n}e^{-j\psi_i} -j \frac{1}{2} \sum_{i=0}^{i=I}\sum_{n=0}^{N-1}x_{n} C_iW_N^{(k-c_i)n}e^{j\psi_i} \hspace{4cm} \\ &=\frac{1}{2} \sum_{n=0}^{N-1} x_{n} \left( \sum_{i=0}^{i=I} j C_iW_N^{(k+c_i)n}e^{-j\psi_i} \right) + \frac{1}{2} \sum_{n=0}^{N-1} x_{n} \left( \sum_{i=0}^{i=I} -j C_iW_N^{(k-c_i)n}e^{j\psi_i} \right) \\ &=\frac{1}{2} \sum_{n=0}^{N-1}x_{n}* \left( \sum_{i=0}^{i=I} C_i*j(W_N^{c_in+N\psi _i /2\pi}-W_N^{-(c_i n+N\psi _i/2\pi)}) \right)* W_N^{kn} \\ &=\frac{1}{2} \sum_{n=0}^{N-1}x_{n}* \left( \sum_{i=0}^{i=I} C_i*j(-2jsin(\frac{2\pi c_i n}{N} + \psi_i)) \right)* W_N^{kn} \\ &=\sum_{n=0}^{N-1}x_{n}* \left( \sum_{i=0}^{i=I} C_i*sin(\frac{2\pi c_i n}{N} + \psi_i) \right)* W_N^{kn} \\ \end{aligned} DFT(Hilbert(xn))=j21i=0i=In=0N1xnCiWN(k+ci)nejψij21i=0i=In=0N1xnCiWN(kci)nejψi=21n=0N1xn(i=0i=IjCiWN(k+ci)nejψi)+21n=0N1xn(i=0i=IjCiWN(kci)nejψi)=21n=0N1xn(i=0i=ICij(WNcin+Nψi/2πWN(cin+Nψi/2π)))WNkn=21n=0N1xn(i=0i=ICij(2jsin(N2πcin+ψi)))WNkn=n=0N1xn(i=0i=ICisin(N2πcin+ψi))WNkn
也就是说被多个余弦信号合成的信号调制后的信号 x n ′ x'_n xn,做希尔伯特变换相当于
H i l b e r t ( x n ′ ) = H i l b e r t ( x n ∗ ∑ i = 0 I C i ∗ c o s ( 2 π c i n N + ψ i ) ) = x n ∗ ∑ i = 0 I C i ∗ s i n ( 2 π c i n N + ψ i ) \begin{aligned} Hilbert(x'_n) &= Hilbert \left( x_{n}*\sum_{i=0}^I C_i * cos(\frac{2\pi c_i n}{N} + \psi_i) \right) \hspace{4cm} \\ &=x_{n}*\sum_{i=0}^I C_i * sin(\frac{2\pi c_i n}{N} + \psi_i) \end{aligned} Hilbert(xn)=Hilbert(xni=0ICicos(N2πcin+ψi))=xni=0ICisin(N2πcin+ψi)

被多个余弦信号合成的信号调制的信号 的复数解析信号

有点绕口,其实这里多个余弦信号合成的信号,相当于一个任意周期性信号,但是它所含有的频率都比原始信号所含的频率都要大。记复数解析信号为 z n z_n zn它可以写成:
z n = x n ′ + j ∗ H i l b e r t ( x n ′ ) = x n ∗ ∑ i = 0 I C i ∗ c o s ( 2 π c i n N + ψ i ) + j ∗ x n ∗ ∑ i = 0 I C i ∗ s i n ( 2 π c i n N + ψ i ) = x n ∗ ( ∑ i = 0 I C i ∗ c o s ( 2 π c i n N + ψ i ) + ∑ i = 0 I j ∗ C i ∗ s i n ( 2 π c i n N + ψ i ) ) = x n ∗ ( ∑ i = 0 I C i ( c o s ( 2 π c i n N + ψ i ) + j s i n ( 2 π c i n N + ψ i ) ) ) = x n ∑ i = 0 I C i e j ( 2 π c i n N + ψ i ) = x n ∑ i = 0 I C i e j ψ i e j 2 π c i n N \begin{aligned} z_n &= x'_n + j*Hilbert(x'_n) \hspace{4cm} \\ &=x_{n}*\sum_{i=0}^I C_i * cos(\frac{2\pi c_i n}{N} + \psi_i) + j*x_{n}*\sum_{i=0}^I C_i * sin(\frac{2\pi c_i n}{N} + \psi_i) \\ &=x_{n}* \left( \sum_{i=0}^I C_i * cos(\frac{2\pi c_i n}{N} + \psi_i) +\sum_{i=0}^I j*C_i * sin(\frac{2\pi c_i n}{N} + \psi_i) \right) \\ &=x_{n}* \left( \sum_{i=0}^I C_i ( cos(\frac{2\pi c_i n}{N} + \psi_i) + jsin(\frac{2\pi c_i n}{N} + \psi_i)) \right) \\ &=x_{n} \sum_{i=0}^I C_i e^{j(\frac{2\pi c_i n}{N} +\psi _i)} \\ &=x_{n} \sum_{i=0}^I C_i e^{j\psi _i} e^{j\frac{2\pi c_i n}{N}} \\ \end{aligned} zn=xn+jHilbert(xn)=xni=0ICicos(N2πcin+ψi)+jxni=0ICisin(N2πcin+ψi)=xn(i=0ICicos(N2πcin+ψi)+i=0IjCisin(N2πcin+ψi))=xn(i=0ICi(cos(N2πcin+ψi)+jsin(N2πcin+ψi)))=xni=0ICiej(N2πcin+ψi)=xni=0ICiejψiejN2πcin

本来想从这个复数解析信号中获得关于包络、瞬时相位的特性,但看看要从 x n ∑ i = 0 I C i e j ψ i e j 2 π c i n N x_{n} \sum_{i=0}^I C_i e^{j\psi _i} e^{j\frac{2\pi c_i n}{N}} xni=0ICiejψiejN2πcin 获取到 x n x_n xn,本人想不到办法。

希尔伯特变换的包络、瞬时相位、瞬时频率概念

由于上面把载波信号假设为多个余弦信号相加而成的,无法得到原始信号,估计是自己的方向错了,于是查了一下资料,找到了一个能看得明白的:
在这里插入图片描述
从离散数据角度,载波信号设置为 c o s ( 2 π c n N + ψ n ) cos(\frac{2\pi c n}{N} + \psi_n) cos(N2πcn+ψn),当然,c大于原始信号 x n x_n xn所含最大频率m又小于 N 2 \frac{N}{2} 2N,而 ψ n \psi _n ψn个数与原始信号个数相等,要求 x n ∗ c o s ψ n x_n*cos\psi_n xncosψn x n ∗ s i n ψ n x_n*sin\psi_n xnsinψn所包含的最大频率也小于c。
再来看看调制后的信号 x n ′ = x n ∗ c o s ( 2 π c n N + ψ n ) x'_n = x_n * cos(\frac{2\pi c n}{N} + \psi_n) xn=xncos(N2πcn+ψn)的傅里叶变换:
X ′ ( k ) = ∑ n = 0 N − 1 x n ∗ W N ( c n + N ψ n / 2 π ) + W N − ( c n + N ψ n / 2 π ) 2 ∗ W N k n = 1 2 ∑ n = 0 N − 1 x n ∗ ( W N c n W ψ n / 2 π + W N − c n W − ψ n / 2 π ) ∗ W N k n = 1 2 ∑ n = 0 N − 1 x n W N ( k + c ) n W ψ n / 2 π + 1 2 ∑ n = 0 N − 1 x n W N ( k − c ) n W − ψ n / 2 π = 1 2 ∑ n = 0 N − 1 x n e − j ψ n W N ( k + c ) n + 1 2 ∑ n = 0 N − 1 x n e j ψ n W N ( k − c ) n = 1 2 ∑ n = 0 N − 1 x n ( c o s ψ n − j s i n ψ n ) W N ( k + c ) n + 1 2 ∑ n = 0 N − 1 x n ( c o s ψ n + j s i n ψ n ) W N ( k − c ) n = 1 2 ( ∑ n = 0 N − 1 x n c o s ψ n W N ( k + c ) n − j ∑ n = 0 N − 1 x n s i n ψ n W N ( k + c ) n ) + 1 2 ( ∑ n = 0 N − 1 c o s ψ n x n W N ( k − c ) n + j ∑ n = 0 N − 1 x n s i n ψ n W N ( k − c ) n ) \begin{aligned} X'(k) &= \sum_{n=0}^{N-1}x_{n}*\frac{W_N^{(cn+N\psi _n/2\pi)}+W_N^{-(cn+N\psi _n/2\pi)}}{2}*W_N^{kn} \hspace{4cm} \\ &=\frac{1}{2} \sum_{n=0}^{N-1}x_{n}*(W_N^{cn}W^{\psi _n/2\pi}+W_N^{-cn}W^{-\psi _n/2\pi})*W_N^{kn} \\ &=\frac{1}{2} \sum_{n=0}^{N-1}x_{n}W_N^{(k+c)n}W^{\psi _n/2\pi} + \frac{1}{2} \sum_{n=0}^{N-1}x_{n}W_N^{(k-c)n}W^{-\psi _n/2\pi} \\ &=\frac{1}{2} \sum_{n=0}^{N-1}x_{n}e^{-j\psi _n}W_N^{(k+c)n} + \frac{1}{2} \sum_{n=0}^{N-1}x_{n}e^{j\psi _n} W_N^{(k-c)n} \\ &=\frac{1}{2} \sum_{n=0}^{N-1}x_{n} (cos\psi_n-jsin\psi_n) W_N^{(k+c)n} + \frac{1}{2} \sum_{n=0}^{N-1}x_{n} (cos\psi _n+jsin\psi _n) W_N^{(k-c)n}\\ &=\frac{1}{2}\left(\sum_{n=0}^{N-1} x_{n} cos\psi_n W_N^{(k+c)n} - j \sum_{n=0}^{N-1}x_{n} sin\psi_n W_N^{(k+c)n} \right) +\frac{1}{2}\left(\sum_{n=0}^{N-1} cos \psi_n x_{n}W_N^{(k-c)n} + j\sum_{n=0}^{N-1}x_{n} sin\psi_n W_N^{(k-c)n} \right) \\ \end{aligned} X(k)=n=0N1xn2WN(cn+Nψn/2π)+WN(cn+Nψn/2π)WNkn=21n=0N1xn(WNcnWψn/2π+WNcnWψn/2π)WNkn=21n=0N1xnWN(k+c)nWψn/2π+21n=0N1xnWN(kc)nWψn/2π=21n=0N1xnejψnWN(k+c)n+21n=0N1xnejψnWN(kc)n=21n=0N1xn(cosψnjsinψn)WN(k+c)n+21n=0N1xn(cosψn+jsinψn)WN(kc)n=21(n=0N1xncosψnWN(k+c)njn=0N1xnsinψnWN(k+c)n)+21(n=0N1cosψnxnWN(kc)n+jn=0N1xnsinψnWN(kc)n)

k+c 和 k-c的特性还是没有改变,所以 1 2 ∑ n = 0 N − 1 x n W − ψ n / 2 π W N ( k − c ) n \frac{1}{2} \sum_{n=0}^{N-1}x_{n}W^{-\psi _n/2\pi}W_N^{(k-c)n} 21n=0N1xnWψn/2πWN(kc)n代表着正半频率频谱, 1 2 ∑ n = 0 N − 1 x n W ψ n / 2 π W N ( k + c ) n \frac{1}{2} \sum_{n=0}^{N-1}x_{n}W^{\psi _n/2\pi}W_N^{(k+c)n} 21n=0N1xnWψn/2πWN(k+c)n代表着负半频率频谱(或者说是大于 N 2 \frac{N}{2} 2N倍频的频谱)

对调制后的信号做希尔伯特变换,再做傅里叶变换,可写成下面这样:

D F T ( H i l b e r t ( x n ′ ) ) = j 1 2 ∑ n = 0 N − 1 x n W N ( k + c ) n W ψ n / 2 π − j 1 2 ∑ n = 0 N − 1 x n W N ( k − c ) n W − ψ n / 2 π = 1 2 ∑ n = 0 N − 1 x n ∗ ( j W N c n W ψ n / 2 π − j W N − c n W − ψ n / 2 π ) ∗ W N k n = ∑ n = 0 N − 1 x n ∗ j W N ( c n + N ψ n / 2 π ) − j W N − ( c n + N ψ n / 2 π ) 2 ∗ W N k n = ∑ n = 0 N − 1 x n ∗ j ( − 2 j s i n ( 2 π c n N + ψ n ) ) 2 ∗ W N k n = ∑ n = 0 N − 1 x n ∗ s i n ( 2 π c n N + ψ n ) ∗ W N k n \begin{aligned} DFT(Hilbert(x'_n)) &=j\frac{1}{2} \sum_{n=0}^{N-1}x_{n}W_N^{(k+c)n}W^{\psi _n/2\pi} -j \frac{1}{2} \sum_{n=0}^{N-1}x_{n}W_N^{(k-c)n}W^{-\psi _n/2\pi} \hspace{4cm} \\ &=\frac{1}{2} \sum_{n=0}^{N-1}x_{n}*(jW_N^{cn}W^{\psi _n/2\pi}-jW_N^{-cn}W^{-\psi _n/2\pi})*W_N^{kn} \\ &=\sum_{n=0}^{N-1}x_{n}*\frac{jW_N^{(cn+N\psi _n/2\pi)}-jW_N^{-(cn+N\psi _n/2\pi)}}{2}*W_N^{kn} \\ &=\sum_{n=0}^{N-1}x_{n}*\frac{j(-2jsin(\frac{2\pi c n}{N} + \psi_n))}{2}*W_N^{kn} \\ &=\sum_{n=0}^{N-1}x_{n}*sin(\frac{2\pi c n}{N} + \psi_n)*W_N^{kn} \\ \end{aligned} DFT(Hilbert(xn))=j21n=0N1xnWN(k+c)nWψn/2πj21n=0N1xnWN(kc)nWψn/2π=21n=0N1xn(jWNcnWψn/2πjWNcnWψn/2π)WNkn=n=0N1xn2jWN(cn+Nψn/2π)jWN(cn+Nψn/2π)WNkn=n=0N1xn2j(2jsin(N2πcn+ψn))WNkn=n=0N1xnsin(N2πcn+ψn)WNkn

也就是说 H i l b e r t ( x n ′ ) = x n ∗ s i n ( 2 π c n N + ψ n ) = x n ∗ c o s ( 2 π c n N + ψ n − π 2 ) Hilbert(x'_n) = x_{n}*sin(\frac{2\pi c n}{N} + \psi_n) = x_{n}*cos(\frac{2\pi c n}{N} + \psi_n - \frac{\pi}{2}) Hilbert(xn)=xnsin(N2πcn+ψn)=xncos(N2πcn+ψn2π)
调制后的信号 的复数解析信号 z n z_n zn为:
z n = x n ′ + j H i l b e r t ( x n ′ ) = x n c o s ( 2 π c n N + ψ n ) + j x n s i n ( 2 π c n N + ψ n ) = x n e j ( 2 π c n N + ψ n ) \begin{aligned} z_n &=x'_n + jHilbert(x'_n)\hspace{4cm} \\ &= x_n cos(\frac{2\pi c n}{N} + \psi_n) + j x_n sin(\frac{2\pi c n}{N} + \psi_n)\\ &=x_n e^{j(\frac{2\pi c n}{N} + \psi_n)} \\ \end{aligned} zn=xn+jHilbert(xn)=xncos(N2πcn+ψn)+jxnsin(N2πcn+ψn)=xnej(N2πcn+ψn)
可以看到,如果 x n x_n xn是正数,对 z n z_n zn求模就可以得到 x n x_n xn了。

如果说TD-OCT接收到的信号通过计算复数解析信号的模可以获得纵深信息那么
估计TD-OCT接收到的信号也是 x n ∗ c o s ( 2 π c n N + ψ n ) x_n * cos(\frac{2\pi c n}{N} + \psi_n) xncos(N2πcn+ψn)形式的。

希尔伯特-黄变换

再深入查了一下 希尔伯特-黄变换:
找到两篇文章描述得十分易懂:
希尔伯特-黄变换(HHT)的前世今生
这篇文章能让你明白经验模态分解(EMD)

网上找的大牛实现的经验模态分解库 PyEMD (不同于pypi上的)
安装方法如下
Installation
Recommended

Simply download this directory either directly from GitHub, or using command line:

$ git clone https://github.com/laszukdawid/PyEMD

Then go into the downloaded project and run from command line:
(使用这种安装方法时 不要把解压后的文件夹放置到python/lib/site-packages)
$ python setup.py install

PyPi

Packaged obtained from PyPi is/will be slightly behind this project, so some features might not be the same. However, it seems to be the easiest/nicest way of installing any Python packages, so why not this one?

$ pip install EMD-signal

测试这个PyEMD库的博客

# 导入工具库
import numpy as np
from PyEMD import EMD, Visualisation

# 构建信号
t = np.arange(0,1, 0.01)
S = 2*np.sin(2*np.pi*15*t) +4*np.sin(2*np.pi*10*t)*np.sin(2*np.pi*t*0.1)+np.sin(2*np.pi*5*t)

# 提取imfs和剩余信号res
emd = EMD()
emd.emd(S)
imfs, res = emd.get_imfs_and_residue() #imf分解 和残留
# 绘制 IMF
vis = Visualisation()
vis.plot_imfs(imfs=imfs, residue=res, t=t, include_residue=True)

在这里插入图片描述
上面模拟的TD-OCT 信号进行HHT变换,求IMF1的包络:

from PyEMD import EMD,Visualisation

#orgsignal = data
orgsignal = amp
#print(np.mean(orgsignal))
submean = np.array(orgsignal,dtype=np.float) - np.mean(orgsignal)
plt.plot(submean)

emd = EMD()
emd.emd(submean)

imfs, res = emd.get_imfs_and_residue() #imf分解 和残留
# 绘制 IMF
vis = Visualisation()
vis.plot_imfs(imfs=imfs, residue=res, include_residue=True)

imf1fft = np.fft.fft(imfs[0])
imf1fft_cnt = len(imf1fft)//2
for i in range(imf1fft_cnt):
    imf1fft[i] = 0-1j * imf1fft[i]
    imf1fft[i+hcnt] = 0+1j * imf1fft[i+hcnt]

imf1fft[0] = 0+0j
imf1fft[imf1fft_cnt] = 0+0j
imf1hilbert = np.fft.ifft(imf1fft)
imf1_z = imfs[0] + 0+1j * imf1hilbert.real

plt.figure()
plt.plot(imfs[0])
plt.plot(imf1hilbert.real)
plt.plot(np.abs(imf1_z))

plt.figure()
plt.plot(np.abs(imf1_z))

在这里插入图片描述
在这里插入图片描述
在这里插入图片描述
在这里插入图片描述
看起来比直接做希尔伯特变换求包络再进行滤波 的效果要好一点。
但是实际的TD-OCT信号不知道是不是效果好一点。

评论 6
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值