numpy 傅立叶得到幅值和频率

做个备份，对 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('
')
    s = s.strip('
')

    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()