加速MATLAB中的傅里叶级数for循环

Question

加速MATLAB中的傅里叶级数for循环

performancematlabfor-loopsumvectorization

3

我知道这个话题已经被多次讨论，提前道歉。我只是无法解决一个效率低下的for循环操作，希望你能帮助我。

我正在实现一个单一的for循环，用于求和时域周期数据。我使用之前从FFT获取的傅里叶系数。我的问题是，我需要添加额外的时域数据到周期信息中，所以我不能简单地使用逆FFT操作。我基本上正在编辑标准离散傅里叶变换（DFT）和，即：

for k = 1:L
    x = x + F(k+1)*exp(1j*(omega(k+1)*t));
end

在这个系列中，F是我的傅里叶系数，ω是频率向量，t是时间向量。我必须修改总和，使其看起来像这样：

for k = 1:L
    x = x + F(k+1)*exp(1j*(omega(k+1)*t)+1j*k*xt);
end

唯一的变化是我现在包括了一个向量，它是一个时间域函数(例如正弦波)，我称之为xt。问题是时间域信息必须有极其细微的分辨率(10秒录音大约长度为5e6)。我不能缩短时间长度，因为我需要高频率分辨率。这导致我的机器上单个函数评估大约需要7小时的时间(我不得不承认，这不是最好的情况)。我需要在优化设置中评估函数，因此7小时的函数评估时间是不可行的。

我尝试过将操作向量化，但矩阵变得太大，我的计算机无法处理，并且MATLAB已经更新以更有效地处理for循环。我尝试编写自己的快速傅里叶变换版本，但由于我在每个步骤中编辑频率信息，某些Cooley和Tukey算法所需的假设就会破裂。有谁知道如何将上述求和重写为更有效的格式吗？我已经预先分配了向量。

- Chris Church

1

好的，这很困难。我不知道我是否能够改进代码。困难的问题是，你真的需要所有这些系数吗？对于许多应用程序来说，丢弃数据的某些方面是很常见的。 - Sanjay Manohar

2

你把求和变成了递归……这是你想要做的吗？

for k = 1:L     x = x + F(k+1)*exp(1j*(omega(k+1)*t)+1j*k*xt);  ## 每次迭代都会更新x，然后在指数中使用x end

- gariepy

这正是我担心的事情。我试图捕捉的现象在整体上非常微小（具体来说是压力信号中的涡轮叶片振动）。叶片振动是多普勒频移的，因此我需要至少32 kHz的频率分辨率以防止混叠。在其当前形式下，我有一个大约52 kHz的上限频率，因此将其降低到32 kHz是一种选择，但我仍然希望尽可能地保留信号。 - Chris Church

是的，我确实打算在每次迭代中添加一个新组件。我尝试根据原始FFT算法以递归格式编写它，但失败了。许多添加部分是重复的，这就是为什么我认为可能改进其当前形式的原因。术语“xt”在每个总和中保持不变。该方程采用标准指数傅里叶级数格式。我的附加组件改变了总和中每个频率分量的相位。 - Chris Church

1

如果你对非常高频的组件感兴趣，那么你是否需要低频呢？换句话说，最慢的fft组件描述了事物如何在非常缓慢的时间尺度上变化。根据频谱F和调整xt所代表的内容，你可能可以忽略其中很多。 - Sanjay Manohar

显示剩余3条评论

2个回答

1

我认为这可能是一种矢量化方法来实现你想要的，使用 bsxfun。

k = 1:L;
x = x + sum(bsxfun(@times, F(k+1), exp(1j*(bsxfun(@times, omega(k+1), permute(t, [2 1])) + bsxfun(@times, k, permute(xt, [2 1]))))), 2).';

注意: 对于这个解决方案，F、omega、k、t和xt都应该是列向量，即1x?。

说明

首先计算指数exp(1j*(omega(k+1)*t)+1j*k*xt)。

e = exp(1j*(bsxfun(@times, omega(k+1), permute(t, [2 1])) + bsxfun(@times, k, permute(xt, [2 1]))))

该程序生成一个 MxL 的矩阵，其中 M 是 t 和 xt 的长度。接下来，我们将这个矩阵的每一列与维度为 1xL 的 F 相乘。

Fe = bsxfun(@times, F(k+1), e)

这是另一个MxL矩阵。最后，我们沿着第二个维度（L）求和，得到一个Mx1向量。

x = x + sum(Fe, 2).'

其中'.'用于将Mx1向量转置为1xM向量。

- IKavanagh

我已经尝试实现您的解决方案。目前我遇到了“两个输入数组的非单例维度必须相互匹配”的错误，这意味着我可能只需要调整数组的维度就可以让它正常工作。如果成功了，我会发布性能改进的结果。也许可以通过在向量化操作之前将数组强制转换为所需的形状来调整解决方案。 - Chris Church

1

@ChrisChurch，我刚刚完成了对我使用的维度的解释。 - IKavanagh

你的解决方案很优雅，毫无疑问它会起作用。然而，当我按照你的建议实施时，我的机器由于工作内存已满而崩溃了。为了应对我可用的8GB RAM，这个解决方案可能需要被“分块”。明天我可能会有一个更大的机器来尝试你的解决方案。我会发布结果的。 - Chris Church

1

@ChrisChurch 谢谢。对于“块状分割”，最好的方法可能是在t和xt上进行。如果您分解了k，则需要存储更多中间数据并执行更多的“sum”操作。祝你好运。希望你成功。 - IKavanagh

非常感谢您以完全向量化的格式提供解决方案！我之前尝试过，但不知道bsxfun，这使得在所需操作下实现严格的矩阵代数解决方案相当困难。如果我能够在更大的机器上直接运行它，我将提供有关结果的信息！我还会给出有关分块和思考的提示。这让我想到其他并行解决方案形式也可能有所帮助（就分块而言）。非常感谢。 - Chris Church

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- user2271770 · Accepted Answer

所以，问题是：

循环速度较慢，但可以解决内存不足的问题；
向量化速度快，但会耗尽你的内存；

~~我无法发布代码，因为你没有提供有关您的 F、omega、t 和 xt 的确切大小和值范围的详细信息，~~但如何使用混合解决方案：将问题分成可以进行矢量化的频率间隔，并对整个频率间隔（块）累积部分解决方案，而不是每个频率。类似于：

    %%
    clear all;  % [!] BEWARE BEFORE EXECUTING, IT DELETES STUFF!!!

    N = 1016419;
    R = 524288;

    t     = linspace(1,10,N);
    F     = 400*rand(1,R) + 200i*rand(1,R);
    omega = 1./(1:R);
    xt    = sin(0.23*t);
    x     = zeros(size(xt));

    U  = 10;
    % ck1 = repmat(1i*t,  U, 1); % trying to be extra clever is
    % ck2 = repmat(1i*xt, U, 1); % not always good for the health
    ck1 = 1i*t;
    ck2 = 1i*xt;
    u0  = 1;
    tic;
    while u0 < R
            u1 = min(u0+U-1,R);

            % x  = x + sum(bsxfun( ...
            %     @times, ...
            %     F(u0+1:u1+1).', ...
            %     exp(omega(u0+1:u1+1).' * ck1 ...
            %       +      (u0+1:u1+1).' * ck2 ...
            %     ) ...
            % ));
            x  = x + sum(bsxfun( ...
                @times, ...
                F(u0+1:u1+1).', ...
                exp(bsxfun(@times, ck1, omega(u0+1:u1+1).') ...
                  + bsxfun(@times, ck2,      (u0+1:u1+1).') ...
                ) ...
            ));
            u0 = u1 + 1;

            fprintf('%d iterations: %.0f sec.\n', u1, toc);
    end;

< p > U 的值应该足够小，以便向量化版本的 xu 贡献适合内存，并且足够大，使向量化有意义。因此，您需要进行一些实验。