# Convolution of two large signals in MATLAB

I need help regarding convolution. I am trying to convolve a signal which is $15\textrm{ ms}$ long, and the sampling frequency is $10\textrm{ GHz}$. I need to convolve it with a rectangular pulse, which has $15\textrm{$\ \mu$s}$ pulse width. When I am trying to use conv function in MATLAB, my PC goes busy state, and it does not give me any answer. Can anybody help me with this problem, please?

• wait, you're really convolving with a rectangle? You're aware that this is just a moving average? Mar 11 '17 at 10:59

If your shorther pulse signal is really a rectengular waveform then please look for the other answer, but otherwise for a general waveform pulse the following code snipped (excerpted from Maximillian's previous post) shows the actual results (on my laptop with MATLAB R2015) of timings for a frequency domain vs time domain implementation of the convolution operation in simplest terms.

Note that I have reduced the sampling rate Fs from 10 Ghz to 1Ghz due to memory reasons and also I have slightly adjusted the sequence lengths to result in a power of $2$ FFT length so that FFT could run at full speed. (this is slightly unfair for a practical course but here I would like to see the maximum attainable speed gain, hence it would be ill posed to use inefficient FFT lengths for such a purpose)

Fs = 10E8; % reduced for MEMORY reasons
T = 17e-3; % slightly adjusted
t = 0:1/Fs:T;
t = t(1:(2^24-2^14+1)); % its length is adjusted to create a 2^N factoring later.
y = sin(2*pi*1000*t); % The input signal
L1 = length(y);

w_rect = 17e-6;  % width of rect
R = ones(1,Fs*w_rect);
R = pulse(1:2^14);
L2 = length(R);

tic
y2 = real(ifft(fft(y,L1+L2-1).*fft(R,L1+L2-1)));
tdft = toc

tic
y3 = conv(y,pulse);
tconv = toc

eff = tconv/tdft


The result is the following: tdft = 1.4769 s, tconv = 42.14 s and hence speed gain = eff = 28.5x.

On the other hand, for slightly different but otherwise unfactorable FFT lengths, the speed gains decreased down to 10x as understandable.

This is exactly the situation where you should seriously consider implementing your convolution in the frequency domain (using efficient FFTs), rather than in the time domain!

The condition is triggered by the fact that the shorter pulse signal (with about 150k samples) is already too long for a direct time domain convolution.

Lets make a brief comparison between the two approaches for the signals of lengths about $L_1 = 150 \times 10^6$ samples and $L_2 = 150 \times 10^3$ samples:

In the first case a direct time domain implementation would approximately require $L_1 \times L_2$ many real MACS which is about $(150 \times 10^6) \times (150 \times 10^3) = 2.25 \times 10^{13}$ real MACS.

Now compare this to an FFT based frequency domain implementation (according to the convolution in time domain equals multiplication in frequency domain theorem of DFT). Again using rough numbers there will be three steps: 1- Convert both signals to frequency domain of the same (and sufficiently long) lengths (which is actually $L_1+L_2-1$ )using efficient FFT, 2-Multpily the obtained DFTs and get the output DFT, 3-Convert the result back to the time domain using the inverse FFT.

The first step requires about $2 \times (150 \times 10^6) \times \log_2(150*10^6) \times 4$ real MACS which is about $3.25 \times 10^{10}$ real MACs

The second step would require roughly $4 \times (150 \times 10^6) = 6 \times 10^8$ real MACs

And the final, third, step would require about $1.6 \times 10^{10}$ real MACS.

Hence the totality of the FFT based frequency domain convolution implementation would require roughly: $(3.2 \times 10^{10}) + (6 \times 10^8) + (1.6 \times 10^{10}) = 5 \times 10^{10}$ real MACS.

Now compare this to the time domain implementation which required $2.25 \times 10^{13}$ real MACS. This could potentially result in about 500x speed improvement assuming a basic implementation. So for example if the convolution would take 5 hours in direct time domain implementation then it should take about 36 seconds in DFT based case!

So you better follow the second approach. That being said, however, now current available multicore CPU + GPU based extensively vectorized (parallel) implementations can dramatically reduce the direct time domain computation time of convolutions, however again then the efficiency of FFT will also benefit from this parallelization and therefore you should, nevertheless, see an order of magnitude improvement in speed when using a DFT approach.

PLease also consider that such long sequences can be processed in blocks and that also would account for certain imprevements in memory and speed performances.

EDIT: I've tried to implement the DFT processing in an non-blocked based mode, where the whole signal processed at once in memory but that required too much RAM (about at least 10 GBs for MATLAB alone in my PC which has only 8GBs of RAM for all programs currently so it started paging and went into freezing, (the computer simply stcuk without any response even to ctrl+alt+del, so I had to hard reset it) If you have more memory such as 16 GB or above you may still try and it may probably do ok, as I assumed it would at most require about 10 GBs of RAM (plus any other programs using your remainig RAM). But otherwise you should consider dividing your larger signal (which 150 M samples) into smaller pieces or use some other technique.

• Thanks for the details explanation. I needed to do in the time domain, I reduced my rectangular pulse for the convolution and got the result. But your explanation helped me a lot for the better understanding. Mar 12 '17 at 7:23

With a general convolution, Fat32 is perfectly right. You should try to perform convolution in the frequency domain, using Overlap-Add or Overlap-Save algorithms (see also here). Actually, I am suprised that MATLAB is not using these algorithms when it comes to convolution (I've tried, and the compuation time was much longer that one would have expected using OLA/OLS algorithms).

However, for your special problem, there is even a simpler solution: The convolution with a rect is actually a moving average filter:

$$y[n] = \sum_{k=n-N}^n x[n]$$ where $N$ is the length of the rect.

We can also write this as

$$y[n] = \sum_{k=0}^nx[n] - \sum_{k=0}^{n-N}x[n]=X[n]-X[n-N]$$

with

$$X[n]=\sum_{k=0}^nx[n]$$

is the cumulative sum of the signal.

Hence, you can implement the convolution as simple as:

X = cumsum(x);
y = x;
y(N:end) = y(N:end) - X(1:end-N+1);


which is by far faster than any "fast-convolution".

• Please note that primary reasons for implementing an overlap save and/or overlap add based convolutions is not the computational efficiency but is the ability of piece by piece block based processing of long sequences. It's only DFT method that can (for long enough blocks) provide any computational gains unless parallel methods are implemented. And I'm also surprized that as of 2017 Matlab still not using GPU accelerated convolutions. Possibly I dont yet know how to invoke it. Mar 11 '17 at 13:13
• Well, I also used an outdated version of Matlab, but still. Regarding OLA/OLS to my understanding it is also used for computational efficiency, once you have to convolve a long with a short sequence: You can use DFTs that are of (roughly) twice the length of the short sequence, and do block-wise processing. At least regarding memory consumption this is more efficient that doing one big FFT for each sequence and one inverse FFT at the end Mar 11 '17 at 13:43
• For memory efficieny you are right, that's one of the reasons of block based processing. However for computational efficiency, short blocks of DFT will not provide much gains, indeed it can even be worse if shorter than a threshold. Consider a block length of $N$ with direct convolution it requires $N^2$ real MACs. But with DFT it aprx requires $(2N)\log_2(2N) \times 4 \times 3 + 8N$ real MACS. So you may find the minimum $N$ for which speed gain = $$\frac{N^2}{(2N)\log_2(2N) \times 4 \times 3 + 8N}$$ is below $1$ or above $1$ A more realistic computation can yield better estimates Mar 11 '17 at 14:22
• Thanks a lot. Your explanation helped me to understand the theory and figure out the solution. Mar 12 '17 at 7:19