Seismic Unix 编程：npfao有什么用？

注意：需已正确安装Seismic Unix（以下简称 SU ）。

当我试图将SU中一个代码用Python实现时，遇到这样一句：

ntfft = npfao(2*nt,4*nt);

虽然知道它是在 2nt 与 4nt 之间选择一个值，但如何选择却不知道，故通过SU编程一探究竟。

SU编程：不接收参数

在任意位置新建一个文件夹，创建 find_out_npfao.c 文件：

#include "su.h"

int main(){
	int ntfft, nt=100;
	ntfft = npfao(2*nt, 4*nt);
	printf("nt %d ntfft %d\n", nt, ntfft);
	return(CWP_Exit());
}

再创建一个 Makefile 文件：

include $(CWPROOT)/src/Makefile.config

D = $L/libcwp.a $L/libpar.a $L/libsu.a
LFLAGS= $(PRELFLAGS) -L$L -lsu -lpar -lcwp -lm $(POSTLFLAGS)
B = .
PROGS = $B/find_out_npfao
INSTALL	:	$(PROGS)
	@-rm -f INSTALL
	@touch $@
$(PROGS):	$(CTARGET) $D 
	-$(CC) $(CFLAGS) $(@F).c $(LFLAGS) -o $@
	@$(MCHMODLINE)
	@echo $(@F) installed in $B
remake	:
	-rm -f $(PROGS) INSTALL
	$(MAKE)

然后在当前目录打开终端，输入 make 回车，程序编译完成。在终端中输入 ./find_out_npfao 即可运行。

但是更换 nt 就要修改 find_out_npfao.c 文件中 nt 的值再编译一次，有些麻烦，可以将程序改为可接收用户指定参数。

SU编程：可接收参数

修改后的 find_out_npfao.c ：

#include "su.h"
/**************** self documentation ******************/
char *sdoc[] = {
" find_out_npfao - pass								",
" find_out_npfao <stdin >stdout [optional parms]	",
" Required Parameters:								",
" 			None									",
" Optional Parameters:								",
" nt=			number of time samples				",
" Notes:		pass								",
NULL};
/******************* end self doc *********************/
int main(int argc, char **argv){
	int nt, ntfft;
	/* hook up getpar */
	initargs(argc, argv);
	requestdoc(1);
	if (!getparint("nt",&nt)) nt=100;
	ntfft = npfao(2*nt, 4*nt);
	printf("nt %d ntfft %d\n", nt, ntfft);
	return(CWP_Exit());
}

 比之前多了很多行，因为用的是 SU 中接收参数的方法，需要在前面加程序说明即 sdoc ，当然也可以用其他方法实现参数的接收。

运行结果：

npfao浅析

快速傅里叶变换（FFT）要求输入的数量是2的幂，但上面运行结果却不是这样，通过阅读源码，了解到 SU 用的不是普通的 FFT，而是一种“质因数”快速傅里叶变换算法（PFA）。

关于 npfao，使用方法：n=npfao(nmin, nmax)，会在 nmin 和 nmax 之间寻找一个最优的进行傅里叶变换的数量 n，其值仅由集合 {2,3,4,5,7,8,9,11,13,16} 中的互素因子相乘得到，因为 n 不会超过 720720 = 5*7*9*11*13*16 。

至于 Python 实现，我选择 numpy.fft.fft。