매우 시끄러운 데이터를위한 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 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])) 

freck1 목록에는 주어진 프로브의 모든 주파수가 포함되어야합니다. 이 목록을 시각적으로 확인했는데 다른 주파수의 양이 매우 높습니다 (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"목록에 넣어야합니다. 하지만 나는 올바른 주파수를 정확하게 식별 할 수없는 방법에 대해 다소 혼란 스럽습니다. 이론은 하나의 주된 주파수 여야합니다. 어떻게하면 모든 잡음으로부터이 모든 데이터의 "메인"주파수를 정확하게 예측할 수 있습니까?

참고로 본 카맨의 와류 거리를 모델링하고 있으며, 나는 소용돌이 흘림의 빈도를 찾고 있습니다. 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])) 

만 alpha가 기본 파수 같이 계산

freq[i] = freq1[i]*alpha 

어디로 주파수 어레이를 계산하는 것이 유용하다

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

에서가 인덱스를 생성

alpha = 1/Ts 

귀하의 샘플링 기간은 Ts입니다. 제 생각에 freq1은 크기가 조정되지 않았기 때문에 주파수 배열이 너무 높습니다.

다른 시간 단계를 사용하여 데이터를 샘플링하는 경우 numpy.interp (예 :)을 사용하여 균등 한 공간 도메인에서 데이터를 보간해야합니다.

주 주파수를 추정하려면 fft로 변환 된 변수가 더 높은 인덱스를 찾아 해당 인덱스를 freq[i]에 연결하십시오.