Неожиданные результаты БПФ с Python

Я анализирую то, что по сути является формой волны дыхания, построенной в 3 различных формах (данные получены из МРТ, поэтому использовалось несколько периодов эхо-сигнала; см. здесь, если вам нужна краткая информация). Вот пара сегментов двух построенных сигналов для некоторого контекста:

введите здесь описание изображения

Для каждой формы волны я пытаюсь выполнить ДПФ, чтобы обнаружить доминирующую частоту или частоты дыхания.

Моя проблема в том, что когда я рисую ДПФ, которые я выполняю, я получаю одну из двух вещей, в зависимости от библиотеки БПФ, которую я использую. Кроме того, ни один из них не является репрезентативным для того, что я ожидаю. Я понимаю, что данные не всегда выглядят так, как мы хотим, но в моих данных явно присутствуют сигналы, поэтому я ожидаю дискретного Преобразование Фурье для получения частотного пика где-то в разумных пределах. Для справки здесь, я бы ожидал от 80 до 130 Гц.

Мои данные хранятся в фрейме данных pandas, причем данные каждого времени эха находятся в отдельной серии. Я применяю выбранную функцию БПФ (см. код ниже) и получаю разные результаты для каждого из них.

Наука (fftpack)

import pandas as pd
import scipy.fftpack

# temporary copy to maintain data structure
lead_pts_fft_df = lead_pts_df.copy()

# apply a discrete fast Fourier transform to each data series in the data frame
lead_pts_fft_df.magnitude = lead_pts_df.magnitude.apply(scipy.fftpack.fft)
lead_pts_fft_df

Числовой:

import pandas as pd
import numpy as np 

# temporary copy to maintain data structure
lead_pts_fft_df = lead_pts_df.copy()

# apply a discrete fast Fourier transform to each data series in the data frame
lead_pts_fft_df.magnitude = lead_pts_df.magnitude.apply(np.fft.fft)
lead_pts_fft_df

Остальная часть соответствующего кода:

ECHO_TIMES = [0.080, 0.200, 0.400] # milliseconds

f_s = 1 / (0.006) # 0.006 = time between samples
freq = np.linspace(0, 29556, 29556) * (f_s / 29556) # (29556 = length of data)

# generate subplots
fig, axes = plt.subplots(3, 1, figsize=(18, 16))

for idx in range(len(ECHO_TIMES)):
    axes[idx].plot(freq, lead_pts_fft_df.magnitude[idx].real, # real part
                   freq, lead_pts_fft_df.magnitude[idx].imag, # imaginary part

    axes[idx].legend() # apply legend to each set of axes

# show the plot
plt.show()

Пост-ТПФ (SciPy fftpack):

введите здесь описание изображения

Пост-ДПФ (NumPy)

введите здесь описание изображения

Вот ссылка на набор данных (в Dropbox), который использовался для создания эти графики, и здесь есть ссылка на Блокнот Юпитер.

РЕДАКТИРОВАТЬ:

Я воспользовался опубликованным советом и взял величину (абсолютное значение) данных, а также построил график с логарифмической осью Y. Новые результаты опубликованы ниже. Похоже, у меня есть некоторый заворот в моем сигнале. Я использую неправильную шкалу частот? Обновленный код и графики приведены ниже.

# generate subplots
fig, axes = plt.subplots(3, 1, figsize=(18, 16))

for idx in range(len(ECHO_TIMES)):
    axes[idx].plot(freq[1::], np.log(np.abs(lead_pts_fft_df.magnitude[idx][1::])),
                   label=lead_pts_df.index[idx], # apply labels
                   color=plot_colors[idx]) # colors
    axes[idx].legend() # apply legend to each set of axes

# show the plot
plt.show()

введите здесь описание изображения


person questionable_code    schedule 19.07.2018    source источник
comment
На самом деле поймал это как раз перед тем, как вы опубликовали. Исправлено в редакции. Спасибо!   -  person questionable_code    schedule 19.07.2018
comment
Каково среднее значение входного ряда? На графиках показаны значения около 0,99, но на метке вертикальной оси указано нормализованное значение. Являются ли графики репрезентативными для фактических данных, передаваемых функциям БПФ?   -  person Warren Weckesser    schedule 19.07.2018
comment
При преобразовании Фурье вас, как правило, не интересуют реальная и мнимая части (если вам действительно не нужна информация о фазе). Таким образом, вы должны учитывать величину комплексных значений. Это энергия на данной частоте, и это поможет вам определить, какие частоты вносят наибольший вклад в сигнал. Другое дело, что сигналы имеют среднее значение, поэтому нулевая частота, т. е. постоянная, является значимой и доминирует над вкладом других частот. Либо вычтите среднее значение сигнала, либо используйте логарифмическую шкалу Y.   -  person jdowner    schedule 19.07.2018
comment
дополнительный вопрос ... эта ссылка, по которой вы ссылаетесь на mriquestions.com, предоставляет образцы файлов набора данных, аналогичные приведенному выше?   -  person Scott Stensland    schedule 19.07.2018
comment
@ScottStensland, скорее всего, нет; мой набор данных взят из проприетарного МРТ-сканера. Я могу предоставить его вам, если это будет полезно   -  person questionable_code    schedule 19.07.2018
comment
да, если бы вы могли опубликовать этот набор данных или аналогичный, это помогло бы нам/мне решить проблему   -  person Scott Stensland    schedule 21.07.2018
comment
@ScottStensland отредактировано со ссылками на набор данных и блокнот Jupyter. Оба должны быть доступны для загрузки из Dropbox.   -  person questionable_code    schedule 25.07.2018
comment
Вы пробовали построить log(abs(fft))?   -  person Cris Luengo    schedule 26.07.2018
comment
@CrisLuengo Я сделал. Пожалуйста, смотрите редактирование выше.   -  person questionable_code    schedule 02.08.2018
comment
Это выглядит так, как ожидалось. ДПФ сигнала с действительным знаком симметрично. Обычно вы смотрите только на первую половину, которая содержит частоты от 0 до Найквиста (== половина частоты дискретизации). В любом случае, теперь вы можете видеть, что ваш компонент DC намного сильнее, чем любой другой компонент, и он скрывал ваши фактические данные при построении линейного графика. Теперь у вас есть пик в ожидаемом месте, или?   -  person Cris Luengo    schedule 02.08.2018
comment
Возможный дубликат Построение быстрого преобразования Фурье в Python   -  person Cris Luengo    schedule 02.08.2018
comment
Эта ссылка может помочь вам узнать больше о том, как построить свой ДПФ.   -  person Cris Luengo    schedule 02.08.2018


Ответы (1)


Я разобрался со своими проблемами.

Крис Луенго очень помог с эта ссылка, которая помогла мне определить, как правильно масштабировать мои интервалы частот и правильно построить ДПФ.

ДОПОЛНИТЕЛЬНО: я забыл учесть смещение, присутствующее в сигнале. Это не только вызывает проблемы с огромным пиком на 0 Гц в ДПФ, но также отвечает за большую часть шума в преобразованном сигнале. Я использовал scipy.signal.detrend(), чтобы устранить это, и получил очень подходящее ДПФ.

# import DFT and signal packages from SciPy
import scipy.fftpack
import scipy.signal

# temporary copy to maintain data structure; original data frame is NOT changed due to ".copy()"
lead_pts_fft_df = lead_pts_df.copy()

# apply a discrete fast Fourier transform to each data series in the data frame AND detrend the signal
lead_pts_fft_df.magnitude = lead_pts_fft_df.magnitude.apply(scipy.signal.detrend)
lead_pts_fft_df.magnitude = lead_pts_fft_df.magnitude.apply(np.fft.fft)
lead_pts_fft_df

Расположите бины по частоте соответствующим образом:

num_projections = 29556
N = num_projections
T = (6 * 3) / 1000 # 6*3 b/c of the nature of my signal: 1 pt for each waveform collected every third acquisition
xf = np.linspace(0.0, 1.0 / (2.0 * T), num_projections / 2)

Затем постройте:

# generate subplots
fig, axes = plt.subplots(3, 1, figsize=(18, 16))

for idx in range(len(ECHO_TIMES)):
    axes[idx].plot(xf, 2.0/num_projections * np.abs(lead_pts_fft_df.magnitude[idx][:num_projections//2]),
                   label=lead_pts_df.index[idx], # apply labels
                   color=plot_colors[idx]) # colors
    axes[idx].legend() # apply legend to each set of axes

# show the plot
plt.show()

введите здесь описание изображения

person questionable_code    schedule 08.08.2018