【问题标题】:Reconstruct Time Series from FFT frequency and strength data using R使用 R 从 FFT 频率和强度数据重建时间序列
【发布时间】:2019-09-03 13:25:12
【问题描述】:

在对 EEG 测量应用傅立叶变换后,我想以绘图的形式将 FFT 的近似值与原始信号进行比较。我必须将 FFT 中的数据(频率和强度)转换回时间序列。 为了转换原始时间序列,我使用了eegkit packageeegfft method。我得到一个频率和幅度列表来近似原始信号。

这里将 FFT 的两个结果显示为缩短的示例:

# Frequency in Hz
freq <- c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10)  

# Strength in uV
ampl <- c(4.1135352, 5.1272713, 3.2069741, 1.5336438, 2.4301334, 1.0974758, 1.8238327, 0.9637886, 1.1401306, 0.2224472)

是否有可以用来从 FFT 近似的频率和幅度数据重构时间序列的包或方法?


编辑

对于原始信号的重构,是否还需要eegfft方法在结果中返回的相位信息?

# Phase shift in range -pi to pi
phase <- c(0.0000000, 1.1469542, -2.1930702, 2.7361738,1.1597980, 2.6118647, -0.6609641, -2.1508755,1.6584852, -1.2906986)

【问题讨论】:

  • 有 11 个频率,但只有 10 个幅度和相位。我想你的意思不是要包括零频率?
  • @JonSpring 是的,好点,零赫兹的频率没有意义。我会改正的

标签: r fft


【解决方案1】:

我希望这样的事情应该可行。

编辑:我已将phases 设置为默认为零(如果丢失且未传递到data_from_fft)。

freq <- c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10)  
ampl <- c(4.1135352, 5.1272713, 3.2069741, 1.5336438, 2.4301334, 1.0974758, 1.8238327, 0.9637886, 1.1401306, 0.2224472)
phase <- c(0.0000000, 1.1469542, -2.1930702, 2.7361738,1.1597980, 2.6118647, -0.6609641, -2.1508755,1.6584852, -1.2906986)
sampl_freq = 1000

data_from_fft <- function(xmin, xmax, sample_freq, 
                          frequencies, amplitudes, phases = 0) {
  x_vals <- seq(xmin, xmax, length.out = sample_freq * (xmax-xmin))
  y_vals <- x_vals * 0
  for (i in seq_along(x_vals)) {
    # Note, I don't understand why the pi/2 phase adjustment is needed here,
    #   but I couldn't get the right answers out eegfft without it... :-(
    y_vals[i] <- sum(amplitudes * sin(2*pi*frequencies * x_vals[i] + phase + pi/2))
  }
  data.frame(x_vals, y_vals)
}

library(tidyverse)

plot_from_FFT <- data_from_fft(0, 1, sampl_freq, freq, ampl, phase)
ggplot(plot_from_FFT, aes(x_vals, y_vals)) +
  geom_line()

现在,让我们看看我们是否可以使用该输出来重构输入:

eegkit::eegfft(plot_from_FFT$y_vals, lower = 1, upper = 20, Fs = sampl_freq) %>% 
  filter(abs(strength) > 0.1)

   frequency  strength  phase.shift
1          1 4.1158607  0.004451123
2          2 5.1177070  1.154553861
3          3 3.2155744 -2.185185998
4          4 1.5319350  2.739953054
5          5 2.4283426  1.173258629
6          6 1.0813858  2.645126993
7          7 1.8323207 -0.644216053
8          8 0.9598727 -2.138381646
9          9 1.1427380  1.685081744
10        10 0.2312619 -1.265466418

是的!这些非常接近输入。

eegkit::eegfft(plot_from_FFT$y_vals, lower = 1, upper = 20, Fs = sampl_freq) %>% 
      filter(abs(strength) > 0.1) %>%
      left_join(
        tibble(frequency = freq,
               strength_orig = ampl,
               phase_orig   = phase)
      ) %>%
      gather(stat, value, -frequency) %>%
      mutate(category = if_else(stat %>% str_detect("str"), "strength", "phase"),
             version  = if_else(stat %>% str_detect("orig"), "plot inputs", "reconstructed inputs"),) %>%
      ggplot(aes(frequency, value, shape = version, size = version)) + 
      geom_point() +
      scale_x_continuous(breaks = 1:10, minor_breaks = NULL) +
      scale_shape_manual(values = c(16, 21)) +
      scale_size_manual(values = c(1,5)) +
      facet_wrap(~category)

【讨论】:

  • 转换矩阵中的代码对我不起作用,因为缺少相位对象
  • 另外,最后一个带有糟糕代码的代码对我不起作用。我对 %>% 运算符不太熟悉。有一个错误告诉我 filter(abs...) 部分缺少 x 值。尽管有这些小缺陷,重建的代码工作正常,所以一旦最后的情节部分工作正常,我就会接受答案。感谢您的帮助
  • %&gt;% 是“管道”运算符。它最初来自magrittr 包,并在tidyverse 包中广泛使用,以更易读的方式将操作步骤链接在一起。它将结果传递给下一个函数的第一项。所以df %&gt;% filter(b&gt;5) %&gt;% rename(x = b) %&gt;% ggplot() 等价于ggplot(rename(filter(df, x&gt;5), x=b))
  • 了解更多关于您的错误的详细信息有助于准确了解正在发生的事情。 library(tidyverse) 在您运行该行之前是否加载无错误? (这将加载代码中使用的包,包括dplyrtidyrggplot2。)运行eegkit::eegfft(plot_from_FFT$y_vals, lower = 1, upper = 20, Fs = sampl_freq) 是否会给出看起来像我得到的输出,但有更多行?确切的错误文本是什么?
  • 顺便说一句,如果您将正弦替换为余弦,则无需调整 pi/2:cos(2*pi*frequencies * x_vals[i] + phase)。显然:cos(x) = sin(x + pi/2)
猜你喜欢
  • 2015-05-20
  • 2013-02-12
  • 1970-01-01
  • 1970-01-01
  • 2011-05-25
  • 2011-09-06
  • 2018-06-24
  • 2021-05-12
  • 2013-10-13
相关资源
最近更新 更多