做个备份，对 numpy 不熟，每次都找函数找半天。
代码里分几块：
1. 从 argc[1] 的文档中读取数据，并转化为 float。文档中有 2001 行，第一行为头，后面 2000 个是采样数据加换行。
2. wd[] 是我这里用的窗函数。是我们某产品里用的窗，参考下面一行 c：
r[p]=s[p] *(0.3125-0.46875*std::cos(2*PI*p/len)+0.1875*std::cos(4*PI*p/len)-0.03125*std::cos(6*PI*p/len));//FDMS_1 Window
3. 使用 numpy 进行 fft。
4. 找出最大和次最大，然后估算真实最大。
5. 画图。
#!/usr/bin/env python # -*- coding: utf-8 -*- import sys import numpy as np import matplotlib.pyplot as pl points1 = 2000 if len(sys.argv) == 1: print "error, no txt file" exit(1) fp = open(sys.argv[1], ‘r‘) s = fp.readline() samples=[] i = 0 wd=[] while i < 2000: wd.append( 0.3125 - 0.46875 * np.cos ( 2 * np.pi * i / 2000 ) + 0.1875 * np.cos( 4 * np.pi * i / 2000 ) - 0.03125 * np.cos( 6 * np.pi * i / 2000 ) ) s = fp.readline() s = s.strip(‘\n‘) s = s.strip(‘\r‘) if s: fl = float(s)*wd[i] samples.append(fl) # print fl else: break i = i + 1 dft_a = np.fft.fft(samples) dft_a_abs = np.fft.fftshift(abs(dft_a)) * 2 / 2000 idx = np.argmax(dft_a_abs) if dft_a_abs[idx - 1] > dft_a_abs[ idx + 1] : idx = idx - 1 y=dft_a_abs[idx] y2=dft_a_abs[idx+1] beta = (y2-y)/(y2+y) alpha = 3.5*beta; rms=(y+y2) * (3.43612+0.85434*pow(alpha,2) + 0.11871*pow(alpha,4) ) / 2 / 1.4142135623 freq = (idx+alpha+0.5)*1000/2000; print rms, freq axis_a = np.linspace(-points1/2, points1/2-1, num=points1) pl.subplot(121) pl.plot(axis_a,dft_a_abs) pl.axis([95, 105, 0, 0.05]) pl.subplot(122) axis_b = np.linspace(0,2000,2000) pl.plot(axis_b,wd) pl.show()
