How to speed up convolution with a million data points

Question

Runzi Hao 2024-5-3

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2114886-how-to-speed-up-convolution-with-a-million-data-points

编辑： Runzi Hao 2024-5-6

I am currently doing convolution using nested for loops, for 10^6 data points in each for loop. Are there ways to speed up the following code? Thanks in advance!

% nIters = 40;
% n = 1e6; 
% mzL = rand(nIters, n); 
% gg  = rand(1, n);
mzR_temp = zeros(nIters, 1);
for c = 1:n 
    mzR_temp(:) = 0;
    for d = 1:c
        mzR_temp(:) = mzR_temp(:) + gg(c-d+1) * mzL(:,d);
    end
    
    mzR_II(:,c) = mzR_temp;
end

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Matt J 2024-5-3

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2114886-how-to-speed-up-convolution-with-a-million-data-points#answer_1452136

编辑：Matt J 2024-5-3

在 MATLAB Online 中打开

Use conv,

mzR_II=conv(gg,mzL,'same');

or FFTs,

mzR_II=ifft( fft(gg,2*n) .* fft(mzL,2*n) , 'symmetric');
mzR_II=mzR_II(1:n);

4 个评论
显示 2更早的评论隐藏 2更早的评论

Paul 2024-5-3

编辑：Paul 2024-5-3

在 MATLAB Online 中打开

It appears that conv has to be 'full' not 'same', after and then truncated.

n = 5;
rng(100)
mzL = rand(1,n);
gg = rand(1,n-1);
for c = 2:n 
    mzR_temp = 0;
    for d = 1:(c-1)
        mzR_temp = mzR_temp + gg(c-d) * mzL(d);
    end
    
    mzR_II(c) = mzR_temp;
end
mzR_II
mzR_II = 1x5
         0    0.0661    0.3983    0.6871    0.6916
<mw-icon class=""></mw-icon>
<mw-icon class=""></mw-icon>
[0 conv(gg,mzL,'same')] 
ans = 1x5
         0    0.6871    0.6916    0.9559    0.7589
<mw-icon class=""></mw-icon>
<mw-icon class=""></mw-icon>
[0 conv(gg,mzL,'full')]
ans = 1x9
         0    0.0661    0.3983    0.6871    0.6916    0.9559    0.7589    0.1194    0.0006
<mw-icon class=""></mw-icon>
<mw-icon class=""></mw-icon>

fftfilt in the Signal Processing Toolbox can be used for the fft approach.

[0 fftfilt(mzL,gg)]
ans = 1x5
         0    0.0661    0.3983    0.6871    0.6916
<mw-icon class=""></mw-icon>
<mw-icon class=""></mw-icon>

Or filter

[0 filter(mzL,1,gg)]
ans = 1x5
         0    0.0661    0.3983    0.6871    0.6916
<mw-icon class=""></mw-icon>
<mw-icon class=""></mw-icon>

Matt J 2024-5-3

编辑：Matt J 2024-5-3

在 MATLAB Online 中打开

There is absolutely no way the computation should take more than 1 second on any computer made within the last 10 years.

n=1e6;
mzL = rand(1,n);
gg = rand(1,n-1);
tic
mzR_II = [0 fftfilt(mzL,gg)];
toc
Elapsed time is 0.174552 seconds.
tic;
mzR_II=ifft( fft(gg,2*n) .* fft(mzL,2*n) , 'symmetric');
mzR_II=mzR_II(1:n);
toc
Elapsed time is 0.145911 seconds.

Runzi Hao 2024-5-6

Awesome! Thanks for the suggestion of using fftfilt.

It turns out that, on my computer, the convolution ran 43 s for a million data points; and with 40 iterations, the total time was 2600+ s. However, the fftfilt function ran only 8 s for the 40 million data points!