Python быстрое преобразование Фурье для очень шумных данных
У меня есть файл с данными о величине скорости и величине завихренности из симуляции жидкости.
Я хочу выяснить, какова частота этих двух наборов данных.
мой код:
# -*- 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 данные о величине скорости.
Моя цель состоит в том, чтобы найти доминирующую частоту для каждого зонда, а затем собрать все и отобразить их в гистограмме.
Мой вопрос о той части, где говорится, что это происходит, когда происходит волшебство. Я считаю, что БПФ уже работает правильно
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". Но я немного озадачен тем, как я не могу правильно определить правильную частоту. Как и должно быть, в теории должна быть одна основная частота. Как я могу правильно оценить "основную" частоту из всех этих данных из всех шумов
Для справки я моделирую вихревую улицу фон Карманна и ищу частоту появления вихрей. https://en.wikipedia.org/wiki/K%C3%A1rm%C3%A1n_vortex_street
Может кто-нибудь помочь мне, как решить эту проблему?
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],