как сделать фильтр верхних частот? - PullRequest
Новая
       48

как сделать фильтр верхних частот?

0 голосов
/ 24 января 2020

У меня есть 3D-матрица данных об уровне моря (время, у, х), и я нашел спектр мощности, взяв квадрат БПФ, но есть низкие частоты, которые действительно доминируют. Я хочу избавиться от этих низких частот, применив фильтр верхних частот ... как бы я go сделал это? Пример набора данных и структуры / кода приведен ниже:

Это набор данных и создание массивов: 

Yearmin = 2018
Yearmax = 2019
year_len = Yearmax - Yearmin + 1.0 # number of years

direcInput = "filepath"
a = s.Dataset(direcInput+"test.nc", mode='r') 

#creating arrays
lat = a.variables["latitude"][:] 
lon = a.variables["longitude"][:] 
time1 = a.variables["time"][:] #DAYS SINCE JAN 1ST 1950
sla = a.variables["sla"][:,:,:] #t, y, x
time = Yearmin + (year_len * (time1 - np.min(time1)) / ( np.max(time1) - np.min(time1))) 


#detrending and normalizing data 
def standardize(y, detrend = True, normalize = True):
    if detrend == True:
        y = signal.detrend(y, axis=0)
    y = (y - np.mean(y, axis=0))
    if normalize == True:
        y = y / np.std(y, axis=0)
    return y

sla_standard = standardize(sla)

print(sla_standard.shape) = (710, 81, 320)


#fft
fft = np.fft.rfft(sla_standard, axis=0)
spec = np.square(abs(fft))

frequencies = (0, nyquist, df)


#PLOTTING THE FREQUENCIES VS SPECTRUM FOR A FEW DIFFERENT SPATIAL LOCATIONS
plt.plot(frequencies, spec[:, 68,85])
plt.plot(frequencies, spec[:, 23,235])
plt.plot(frequencies, spec[:, 39,178])
plt.plot(frequencies, spec[:, 30,149])
plt.xlim(0,.05)
plt.show()

plot of frequencies and power spectrum

Моя цель - создать фильтр верхних частот ОРИГИНАЛЬНОГО временного ряда (sla_standard), чтобы убрать два действительно больших пика. Какой тип фильтра я должен использовать? Спасибо!

Ответы [ 2 ]

0 голосов
/ 24 января 2020

Различные фильтры: Когда вы заявляете, фильтр действует на частоты вашего сигнала. Существуют различные формы фильтров, и некоторые из них имеют сложные выражения, потому что они должны быть реализованы в режиме реального времени (причинно). Однако в вашем случае вы, похоже, публикуете данные. Вы можете использовать преобразование Фурье, для которого требуются все данные (не причинно-следственные).

Фильтр на выбор: Следовательно, вы можете напрямую выполнить операцию фильтрации в области Фурье, применив маску на ваших частотах. Если вы хотите удалить частоты, я рекомендую вам использовать двоичную маску, состоящую из 0 и 1. Почему? Потому что это самый простой фильтр, о котором вы можете подумать. С научной точки зрения важно заявить, что вы полностью удалили некоторые частоты (скажите это и обоснуйте). Однако сложнее утверждать, что вы позволили некоторым и немного ослабили другие, и что вы произвольно выбрали коэффициент ослабления ...

Python реализация 

signal_fft = np.fft.rfft(sla_standard,axis=0)
mask = np.ones_like(sla_standard)
mask[freq_to_filter,...] = 0.0 # define here the frequencies to filter
filtered_signal = np.fft.irfft(mask*signal_fft,axis=0)
0 голосов
/ 24 января 2020

Используйте .axes.Axes.set_ylim для установки предела по оси Y.

Axes.set_ylim(self, left=None, right=None, emit=True, auto=False, *, ymin=None, ymax=None)

Так что в вашем случае ymin=None, и вы установите ymax, например, на ymax=60000 перед началом построения графика.

Таким образом plt.ylim(ymin=None, ymax=60000).

Извлечение данных не следует делать здесь, потому что это "фальсификация результатов". На самом деле вы хотите увеличить масштаб графика. Человек, который читает таблицу независимо от вас, будет неверно истолковывать данные, если об этом заранее не уведомят. Пики, которые go вне графика - это нормально, потому что все это понимают.

Или:

Непосредственная замена определенных значений в массиве (обр):

arr[arr > ori] = dest

Например, в вашем случае ori=60000 и dest=1

Все значения больше ">", чем 60k, заменяются на 1.

...