numpy 傅立叶得到幅值和频率

时间：2022-11-28 09:09:23

做个备份，对 numpy 不熟，每次都找函数找半天。

代码里分几块：

1. 从 argc[1] 的文档中读取数据，并转化为 float。文档中有 2001 行，第一行为头，后面 2000 个是采样数据加换行。

2. wd[] 是我这里用的窗函数。是我们某产品里用的窗，参考下面一行 c：

r[p]=s[p] *(0.3125-0.46875*std::cos(*PI*p/len)+0.1875*std::cos(*PI*p/len)-0.03125*std::cos(*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()

相关文章