码迷,mamicode.com
首页 > 其他好文 > 详细

numpy 傅立叶得到幅值和频率

时间:2018-01-05 15:16:17      阅读:174      评论:0      收藏:0      [点我收藏+]

标签:bsp   oat   strip   bin   line   ==   fft   nump   max   

做个备份,对 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()

 

numpy 傅立叶得到幅值和频率

标签:bsp   oat   strip   bin   line   ==   fft   nump   max   

原文地址:https://www.cnblogs.com/pied/p/8203693.html

(0)
(0)
   
举报
评论 一句话评论(0
登录后才能评论!
© 2014 mamicode.com 版权所有  联系我们:gaon5@hotmail.com
迷上了代码!