-1

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

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

我的代码：

# -*- 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 containig all the timesteps, cutting of the first 180, because the system has to stabelize.
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 dimensionallity 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, wich all have the data from timestep 180-301
    for l in probeU1:
      b.append(l[i])
      #the freqeuncy was not oscillating around 0, so moved it there by substracting 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('velocitiy 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 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(freq1[val]))

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

作为参考，我正在模拟 Von Karmann 涡街，并且正在寻找涡流脱落的频率。https://en.wikipedia.org/wiki/K%C3%A1rm%C3%A1n_vortex_street

谁能帮我解决这个问题？

4

1 回答 1

0

线

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].

于 2017-05-09T19:18:45.267 回答