Python Frequency filtering с, казалось бы, неправильными частотами

Сценарий ниже фильтрует частоты, обрезая все частоты больше 6. Однако вместо использования, казалось бы, правильной функции rfftfreq, fftfreq используется.

Насколько я понимаю rfftfreq следует использовать вместе с rfft, Почему этот код работает, хотя он использует fftfreq с rfft?

import numpy as np
from scipy.fftpack import rfft, irfft, fftfreq

time   = np.linspace(0,10,2000)
signal = np.cos(5*np.pi*time) + np.cos(7*np.pi*time)

W = fftfreq(signal.size, d=time[1]-time[0])
f_signal = rfft(signal)

# If our original signal time was in seconds, this is now in Hz    
cut_f_signal = f_signal.copy()
cut_f_signal[(W<6)] = 0

cut_signal = irfft(cut_f_signal)

Фон: rfft дает массив, сортирующий моды Фурье с реальным и воображаемым в отдельных записях. Такие как [R(0), R(1), I(1), ... R(N/2),I(N/2)] за R(n) а также I(n) будучи реальной и мнимой частью режима Фурье соответственно. (при условии четного массива записей)

Следовательно, rfftfreq возвращает массив частот, соответствующий этому массиву, например (при условии, что массив четных входов и интервал выборки равен 1):

[0, 1/n, 1/n, 2/n, 2/n ... n/(2n), n/(2n)]

Однако этот код работает с fftfreq где вывод функции

[0, 1/n, ... n/(2n), -n/(2n), ..., -1/n]

Очевидно, что fftfreq должен привести к неправильным результатам при использовании вместе с rfft потому что частоты и бункеры FFT не совпадают.

1 ответ

Решение

Вы неправильно указали частоты в исходном сигнале.

Синусоидальная волна параметризована согласно этому уравнению (из Википедии):

Фактор 2 отсутствует в определении signal = np.cos(5*np.pi*time) + np.cos(7*np.pi*time), Поэтому фактические частоты

5*pi*t = 2*pi*t * f
f = (5*pi*t) / (2*pi*t) = 5 / 2

7*pi*t = 2*pi*t * f
f = (7*pi*t) / (2*pi*t) = 7 / 2

Словом, две частоты - это половина того, что вы думаете. По иронии судьбы, именно поэтому он, кажется, работает с fftfreq вместо rfftfreq, Первый охватывает в два раза частотный диапазон (положительные и отрицательные частоты) и, следовательно, компенсирует отсутствующий фактор 2.

Это правильный код:

signal = np.cos(5 * 2*np.pi * time) + np.cos(7 * 2*np.pi * time)
W = rfftfreq(signal.size, d=time[1]-time[0])
Другие вопросы по тегам