对变形测量信号进行傅里叶变换后，如何确定每个点的频率

变形分析中，常需对类周期信号进行傅里叶变换，以将时域信号转化为频率信号。变换后，如何确定每个点的频率值，相比其振幅值，要复杂一点。

采样频率

采样频率，Fs，也就是单位时间内的采样个数。FFT多用于高频信号的处理中，因此多数帖子中定义采样率为“一秒内采样个数”，可是对于变形监测信号，监测间隔少则一天，多则一月，因此，对于测量信号，建议将采样频率定义为每年采样个数。

计算方法

1. 对变形离散信号S进行FFT变换，方法很多，一个函数而已。
2. 定义采样频率Fs，（每年监测多少期）
3. 用采样频率计算监测信号中周期个数：
   $$Tn=N/Fs$$
   式中，N 为时序 S 的长度 。
4. 确定每个点对应的频率
   $$w=K/Tn$$
   式中，K 为时序S对应的自然数序列，从 0 到 N-1 的自然数序列。

举例

如图所示，ITRS公布的地球极移 X 的10年监测数据。其中，监测间隔为天。
对其进行FFT变换，有：

da1=get_Pxy(da,2000,1,1,2006,12,31)  #获得极移中1997.01.01-2006.12.31的部分
# 对da1x进行傅里叶变换
from scipy.fftpack import fft  #fft表示快速傅里叶变换，ifft表示其逆变换  
N = len(da1x)
fft_y=np.fft.rfft(da1x)  #变换进行FFT
abs_y=np.abs(fft_y)/N               # 取复数的绝对值，即复数的模，获得振幅值,归一化处理
abs_y_half =abs_y[range(int(N/2))]  # 获得单边频谱
# 确定频率。
Fs=365 #将一年作为一个单位时间，该单位时间内采样个数为365个。所以采样率为365
T=N/Fs #用采样率算出段数据中一共有多少个周期
K=np.arange(N)               # 把采样点数的等差数列k除以周期T，就是频率 frq = k/T
freq=K/T                     # 计算每个点的频率值
freq_half=freq[range(int(N/2))]    # 由于对称性，取一半即可

fig3=plt.figure('单边频谱图')
plt.plot(freq_half,abs_y_half,'b-.')
plt.title('单边频谱(归一化)',fontsize=15,color='blue')
plt.xlabel('频率',fontsize=15)
plt.ylabel('振幅',fontsize=15)
plt.show()

为提取有效频率分量，对生成的图片进行局部显示：

fig4=plt.figure('局部显示')
plt.plot(freq_half[1:30],abs_y_half[1:30],'r-^')
plt.title('单边频谱(归一化)',fontsize=15)
plt.xlabel('频率',fontsize=15)
plt.ylabel('振幅',fontsize=15)
plt.show()

上述频率横轴以年为单位，其中两个振幅较大值对应的频率约为1年和0.83年，也即是极移分量中的周年项与钱德勒项。（多数文献提及的半周年分量，不知为何没有显示出来）。
当面对不同形变信号时，需要重新考虑自己定义的“单位时间”，因为要满足萘奎斯特采样定理等等。

这就是最近困扰我的测量信号FFT变换后，每点的频率确定问题。如有错误，欢迎留言指正。