Python 快速傅里叶变换适用于非常嘈杂的数据

Question

我有一个文件，其中包含来自流体模拟的速度幅度数据和涡度幅度数据。

我想知道这两个数据集的频率是多少。

我的代码：

# -*- coding: utf-8 -*-
"""
Spyder Editor

This is a temporary script file.
"""
import re
import math
import matplotlib.pyplot as plt

import numpy as np


probeU1 = []
probeV1 = []
# this creates an array containing all the timesteps, cutting of the first 180, because the system has to stabilize.
number = [ round(x * 0.1, 1) for x in range(180, 301)] 

# this function loops over the different time directories, and reads the velocity file.
for i in range(len(number)):
    filenamepath = "/Refinement/Vorticity4/probes/" +str(number[i]) + "/U"
    data= open(filenamepath,"r")
    temparray = []
    #removes all the formatting around the data
    for line in data:
        if line.startswith('#'):
            continue
    else:
        line = re.sub('[()]', "", line)
        values = line.split()
        #print values[1], values[2]
        xco = values[1::3]
        yco = values[2::3]        
        #here it extracts all the velocity data from all the different probes
        for i in range(len(xco)):
            
            floatx = float(xco[i])
            floaty = float(yco[i])
            temp1 = math.pow(floatx,2)
            temp2 = math.pow(floaty,2)
            #print temp2, temp1
            temp3 = temp1+temp2
            temp4 = math.sqrt(temp3)
            #takes the magnitude of the velocity
            #print temp4
            temparray.append(temp4)
        probeU1.append(temparray)
        
#        
#print probeU1[0]  
#print len(probeU1[0])     
#        
        
# this function loops over the different time directories, and reads the vorticity file.
for i in range(len(number)):
    filenamepath = "/Refinement/Vorticity4/probes/" +str(number[i]) + "/vorticity"
    data= open(filenamepath,"r")
#    print data.read()
    temparray1 = []

    for line in data:
        if line.startswith('#'):
            continue
    else:
        line = re.sub('[()]', "", line)
        values = line.split()
        zco = values[3::3]
        #because the 2 dimensionality the z-component of the vorticity is already the magnitude
        for i in range(len(zco)):  
            abso = float(zco[i])
            add = np.abs(abso)
            
            temparray1.append(add)
    
    probeV1.append(temparray1)
    
#Old code block to display the data and check that it made a wave pattern(which it did) 
##Printing all probe data from 180-300 in one graph(unintelligible)
#for i in range(len(probeU1[1])):
#    B=[]
#    for l in probeU1:
#      B.append(l[i])
##    print 'B=', B
##    print i
#    plt.plot(number,B)
#
#
#plt.ylabel('magnitude of velocity')
#plt.show()  
#
##Printing all probe data from 180-300 in one graph(unintelligible)
#for i in range(len(probeV1[1])):
#    R=[]
#    for l in probeV1:
#      R.append(l[i])
##    print 'R=', R
##    print i
#    plt.plot(number,R)
#
#
#plt.ylabel('magnitude of vorticity')
#plt.show()  
    
#Here is where the magic happens, (i hope)
ans=[]
for i in range(len(probeU1[1])):
    b=[]
   #probeU1 is a nested list, because there are 117 different probes, which all have the data from timestep 180-301
    for l in probeU1:
      b.append(l[i])
      #the frequency was not oscillating around 0, so moved it there by subtracting the mean
      B=b-np.mean(b)
      #here the fft happens
      u    = np.fft.fft(B)
      #This should calculate the frequencies?
      freq = np.fft.fftfreq(len(B), d= (number[1] - number[0]))
      # If im not mistakes this finds the peak frequency for 1 probe and passes it another list
      val = np.argmax(np.abs(u))
      ans.append(np.abs(freq[val]))    
 
      plt.plot(freq, np.abs(u))

#print np.mean(ans)
plt.xlabel('frequency?')
plt.savefig('velocity frequency')
plt.show()

# just duplicate to the one above it
ans1=[]
for i in range(len(probeV1[1])):
    c=[]
    
    for l in probeU1:
      c.append(l[i])
      C=c-np.mean(c)
      y    = np.fft.fft(C)
      freq1 = np.fft.fftfreq(len(C), d= (number[1] - number[0]))
      val = np.argmax(np.abs(y))
      ans1.append(np.abs(freq1[val]))    
 
      plt.plot(freq1, np.abs(y))

#print np.mean(ans1)
plt.ylabel('frequency?')
plt.savefig('vorticity frequency')
plt.show()     



data.close()

我的数据包含 117 个探测器，每个探测器都有自己的 121 个速度幅度数据。

我的目标是找到每个探针的主频率，然后收集所有这些频率并将它们绘制在直方图中。

我的问题是关于它说这就是魔法发生的地方的部分。我相信 fft 已经正常工作了

  y    = np.fft.fft(C)
          freq1 = np.fft.fftfreq(len(C), d= (number[1] - number[0]))

如果我没有记错的话，freq1 列表应该包含给定探头的所有频率。我已经目视检查过此列表，不同频率的数量非常高（20+），因此信号可能非常嘈杂。

# If I'm not mistakes this finds the peak frequency for 1 probe and passes it another list
      val = np.argmax(np.abs(u))
      ans.append(np.abs(freq1[val]))    

理论上，这部分应该从一个探头获取最大信号，然后放入“ans”列表中。但我对如何无法正确识别正确的频率感到有点困惑。因为理论上应该有一个主频率。如何从所有噪声的所有数据中正确估计“主”频率

作为参考，我正在对冯卡曼涡街进行建模，并且正在寻找涡旋脱落的频率。 https://en.wikipedia.org/wiki/K%C3%A1rm%C3%A1n_vortex_street

谁能帮我解决这个问题吗？

Answer 1

线路

freq1 = np.fft.fftfreq(len(C), d= (number[1] - number[0]))

仅生成来自

的索引

freq1 = [0, 1, ...,   len(C)/2-1,     -len(C)/2, ..., -1] / (d*len(C))

这对于计算频率数组很有用

freq[i] = freq1[i]*alpha

其中

alpha

是基本波数，计算公式为

alpha = 1/Ts

是

Ts

您的采样周期。我认为因为

freq1

没有缩放，所以频率数组如此之高。

请注意，如果您使用不同的时间步长对数据进行采样，则需要使用

numpy.interp

（例如）在均匀空间域中对其进行插值。

要估计主频率，只需找到 fft 变换变量较高的索引，并将该索引与

freq[i]

相关联。