1

我正在尝试在 R 中复制下图：（改编自http://link.springer.com/article/10.1007/PL00011669）

DFT 图

该图的基本概念是显示 DFT 的前几个分量，在时域中绘制，然后在时域中显示仅使用这些分量 (X') 相对于原始数据 (X) 的重构波。我想稍微修改上图，使显示的所有线条都覆盖在一个图上。

我一直在尝试用一些以 60 Hz 采样的真实数据来调整这个数字。例如：

## 3 second sample where: time is in seconds and var is the variable of interest
temp = data.frame(time=seq(from=0,to=3,by=1/60),
            var = c(0.054,0.054,0.054,0.072,0.072,0.072,0.072,0.09,0.09,0.108,0.126,0.126,
                  0.126,0.126,0.126,0.144,0.144,0.144,0.144,0.144,0.162,0.162,0.144,0.126,
                  0.126,0.108,0.144,0.162,0.18,0.162,0.126,0.126,0.108,0.108,0.126,0.144,
                  0.162,0.144,0.144,0.144,0.144,0.162,0.162,0.126,0.108,0.09,0.09,0.072,
                  0.054,0.054,0.054,0.036,0.036,0.018,0.018,0.018,0.018,0,0.018,0,
                  0,0,-0.018,0,0,0,-0.018,0,-0.018,-0.018,0,-0.018,
                  -0.018,-0.018,-0.018,-0.036,-0.036,-0.054,-0.054,-0.072,-0.072,-0.072,-0.072,-0.072,
                  -0.09,-0.09,-0.108,-0.126,-0.126,-0.126,-0.144,-0.144,-0.144,-0.162,-0.162,-0.18,
                  -0.162,-0.162,-0.162,-0.162,-0.144,-0.144,-0.144,-0.126,-0.126,-0.108,-0.108,-0.09,
                  -0.072,-0.054,-0.036,-0.018,0,0,0,0,0.018,0.018,0.036,0.054,
                  0.054,0.054,0.054,0.054,0.054,0.054,0.054,0.054,0.054,0.072,0.054,0.072,
                  0.072,0.072,0.072,0.072,0.072,0.054,0.054,0.054,0.036,0.036,0.036,0.036,
                  0.036,0.054,0.054,0.072,0.09,0.072,0.036,0.036,0.018,0.018,0.018,0.018,
                  0.036,0.036,0.036,0.036,0.018,0,-0.018,-0.018,-0.018,-0.018,-0.018,0,
                  -0.018,-0.036,-0.036,-0.018,-0.018,-0.018,-0.036,0,0,-0.018,-0.018,-0.018,-0.018))

##plot the original data
ggplot(temp, aes(x=time, y=var))+geom_line()

我相信我可以用它fft()来最终实现这个目标，但是从输出fft()到我的目标的飞跃有点不清楚。

我意识到这个问题有点类似于：如何计算来自实值输入的 fft() 输出的幅度和相位角？但我对上述特定数据的实际代码更感兴趣。

请注意，我对时间序列分析相对较新，因此您可以将 fft() 的输出放在上下文中提供的任何清晰度，或者您可以推荐的任何可以有效完成此任务的包，将不胜感激。

谢谢

4

2 回答 2

1

matlab是你最好的工具，具体的函数就是fft()。要使用它，首先要确定您的时域数据的几个基本参数：

1，持续时间（T），等于3s。

2、采样间隔T_s，等于1/60 s。

3、频域转数f_s，等于相邻两个傅里叶基的频率差。您可以根据需要定义 f_s。然而，最小可能的 f_s 等于 1/T=0.333 Hz。因此，如果您想要更好的频域转数（更小的 f_s），则需要更长的时域数据。

4、最大频率f_M，根据香农采样理论，等于1/(2T_s)=30。

5、DFT长度N，等于2*f_M/f_s。

然后找出要用于近似数据的四个傅里叶基的特定频率。例如，3、6、9 和 12 Hz。所以 f_s = 3 Hz。那么N=2*f_M/f_s=20。

您的 Matlab 代码如下所示：

var=[0.054,0.054,0.054 ...]; % input all your data points here
f_full=fft(var,20); % Do 20-point fft
f_useful=f_full(2:5); % You are interested with the lowest four frequencies except DC

这里 f_useful 包含四个傅里叶基的四个复系数。要重建 var，请执行以下操作：

% Generate basis functions
dt=0:1/60:3;
df=[3:3:12];
basis1=exp(1j*2*pi*df(1)*dt);
basis2=exp(1j*2*pi*df(2)*dt);
basis3=exp(1j*2*pi*df(3)*dt);
basis4=exp(1j*2*pi*df(4)*dt);

% Reconstruct var
var_recon=basis1*f_useful(1)+...
basis2*f_useful(2)+...
basis3*f_useful(3)+...
basis4*f_useful(4);
var_recon=real(var_recon);

% Plot both curves
figure;
plot(var);
hold on;
plot(var_recon);

将此代码改编为您的论文:)

于 2013-08-21T02:29:57.487 回答
0

从Signal Processing改编我自己的帖子。我认为它仍然与 Python 相关。

我不是这个主题的专家，但有一些有用的例子可以分享。

您保留的傅立叶分量越多，您就越能模仿原始信号。

此示例显示了当您保留 10、20、...最多n 个组件时会发生什么。假设x和y是您的数据向量。

import numpy
from matplotlib import pyplot as plt

n = len(y)
COMPONENTS = [10, 20, n]

for c in COMPONENTS:
    colors = numpy.linspace(start=100, stop=255, num=c)
    for i in range(c):
        Y = numpy.fft.fft(y)
        numpy.put(Y, range(i+1, n), 0.0)
        ifft = numpy.fft.ifft(Y)
        plt.plot(x, ifft, color=plt.cm.Reds(int(colors[i])), alpha=.70)

    plt.title("First {c} fourier components".format(c=c))
    plt.plot(x,y, label="Original dataset", linewidth=2.0)
    plt.grid(linestyle='dashed')
    plt.legend()
    plt.show()

对于本书的数据集，最多保留 4、10 和n 个分量：

前 4 个傅立叶分量 前 10 个傅里叶分量 前 98 个傅立叶分量

对于您的数据集，最多保留 4、10 和n 个分量：

前 4 个傅立叶分量 前 10 个傅里叶分量 前 181 个傅立叶分量

于 2018-05-16T19:37:28.550 回答