pythonrfft不能产生正确的相位信息

Question

所以我试图写一个代码来实现Forman相位校正，以纠正干涉图中与相位相关的不对称。这个程序的前提是将exp（-i*PhaseAngle）的傅里叶变换与干涉图卷积。这项技术从60年代就已经出现了，应该可以用，但我一直有困难。在

像python这样的干涉图正好用来调试相位I。一旦我开始工作，我就生成了一个相位为0.3弧度的干涉图，而我的相位阵列看起来不像是一个常数0.3，我不确定会发生什么。我尝试过numpy arctan、numpy arctan2和numpy angle，结果都很糟糕。我也尝试过修改我的数据，但这没用，所以为了清晰起见，我删除了它。在

你知道为什么我不能从程序中得到正确的阶段信息吗？如果将这些代码插入python中，就会得到一些不稳定的阶段图。在

我对我的处理函数很有信心，代码的真正内容在下面，我调用这些函数来处理我的干涉图。计算相角的基本程序是：

-将干涉图截短到几个点（我用401） -旋转它，使零路径差（ZPD）成为第一个数据点 -将ZPD值乘以0.5（不确定原因，但学术论文说你必须这样做） -取结果数据的实（单侧）fft -取虚数除以实数的反正切，得到相角

所以我按照上面的程序来制作一个相移为0.3弧度的人工产生的干涉图，但是无法用这个程序重现这个相位误差。有什么想法吗？谢谢你们的帮助！在

编辑：这个问题归根结底是如何安排我的信号（它有一个已知的，恒定的，相位误差为0.3）在应用fft之前，我的结果的虚部是一个常数0.3？而且，我应该在python中使用什么类型的fft来实现这个结果？我非常有信心，我成功地产生了一个恒定相位误差为0.3的信号。在

代码：

## Importer
import numpy as np
import scipy.integrate
import matplotlib.pyplot as plt
plt.style.use('ggplot')
from scipy.fftpack import dct,idct,dst,idst
import cmath
from scipy import signal
from scipy.integrate import simps
import csv

## Generate interferogram
xlist = np.arange(-0.03, 0.03,(1/15798.69))
numPhase = 200
print 'generating interferogram...'
def infGen(nu, xval):
    h = 6.626*10**-34 #J*s
    c = 3*10**8 #m/s
    kB = 1.38*10**-23 #J/K
    T = 300.0 # K    
    return 0.5*2*h*c*c*nu*nu*nu*(100.0**3)*1/(np.exp(h*c*nu*100/(kB*T))-1)*np.cos(2*np.pi*nu*xval+0.3)
interferogram = list()
intError = list()
for i in range(len(xlist)):
    xval = xlist[i]
    value, error = scipy.integrate.quad(infGen,500,4000,args = (xval),limit=150)
    interferogram.append(value)
    intError.append(error)

## Processing functions
def truncate(interferogram, centerburst):
    truncated=list()
    truncatedLength=2*numPhase+1
    for i in range(truncatedLength):
        truncated.append(interferogram[centerburst-(truncatedLength/2)+i])
    #Compute interferogram averages to normalize
    truncavg=np.average(truncated)
    #Normalize the interferograms by their averages
    for i in range(len(truncated)):
        truncated[i]=truncated[i]-truncavg
    return truncated
def rotateMertz(interferogram):
    rotated=list()
    for i in range(len(interferogram)-numPhase):
        rotated.append(interferogram[i+numPhase])
    for i in range(numPhase):
        rotated.append(interferogram[i])
    return rotated   
def getphaseangle(fftdata):
    phaseangle=list()
    for i in range(len(fftdata)):
        #phaseangle.append(np.angle(fftdata[i]))
        """if fftdata[i].imag < 0:
            phaseangle[i] = phaseangle[i]+2*np.pi
        if fftdata[i].real>0 and np.abs(fftdata[i].imag)<(10**-10) and fftdata[i].imag<0:
            phaseangle[i] = 0"""
        phaseangle.append(np.arctan(fftdata[i].imag/fftdata[i].real))
    return phaseangle
def findZPD(inf):
    # finds the index of the zpd for a given interferogram
    maxVal=max(np.abs(inf))
    x = [i for i, j in enumerate(inf) if j==maxVal]
    return x[0]

## Process Interferogram
inf = interferogram
phasedataTrunc = truncate(inf, findZPD(inf))
phasedataRot = rotateMertz(phasedataTrunc) # Rotates so ZPD is first value
phasedataRot[0] = phasedataRot[0]*0.5 # Multiply ZPD value by 0.5
phasedataFFT = np.fft.rfft(phasedataRot)
plt.figure('real part')
plt.plot(phasedataFFT.real)
plt.figure('imaginary part')
plt.plot(phasedataFFT.imag)
plt.figure('imaginary/real')
plt.plot(phasedataFFT.imag/phasedataFFT.real)
phaseAngle = getphaseangle(phasedataFFT)
plt.figure('phase angle')
plt.plot(phaseAngle)
plt.show()