Convolution Calculator in MATLAB

Hello friends! Today, we will build a convolution calculator in MATLAB and understand the calculation behind its output. We will first work through two short sequences by hand, then implement the same operation with loops and compare our answer with MATLAB's built-in conv function.

Convolution appears in digital filters, communication systems, image processing and the analysis of linear systems. It combines shifted, weighted contributions from two sequences. Multiplying corresponding samples once is not enough; each output sample contains a sum of products taken at different relative positions.

Our calculator will accept finite, real-valued numerical vectors and their starting sample indices. It will return the full linear convolution, display correctly indexed plots and report the numerical difference between the manual calculation and MATLAB. You can explore related exercises in our MATLAB projects collection.

What Is Discrete-Time Convolution?

Let x[n] represent an input sequence and h[n] represent the impulse response of a linear, time-invariant system. Their convolution gives the zero-state output y[n]. Zero-state means we are calculating the response caused by the input, without a separate contribution from energy initially stored in the system.

y[n] = sum over k of x[k]h[n − k]

For each selected output index n, we multiply samples of x by the corresponding reversed and shifted samples of h, then add the products. Samples outside either finite sequence are treated as zero. Changing n changes which samples overlap, which is why a convolution usually produces more samples than either input.

You can also interpret this calculation without physically reversing a drawing. Each sample x[k] scales a copy of h and shifts that copy by k samples. Adding all those shifted copies produces the output. This second interpretation leads directly to the two nested loops we will use below.

Convolution is useful beyond physical systems. Polynomial multiplication follows the same coefficient calculation, while a moving-average filter convolves a signal with a short sequence of equal weights. The mathematics is shared, but the meaning and units of the numbers depend on the application.

Output length and starting index

If x contains M samples and h contains N samples, their full linear convolution has M + N − 1 samples. This counts every position with possible overlap, including the beginning and ending portions where only a few terms contribute.

If the first sample of x has index n_x and the first sample of h has index n_h, the first output index is n_x + n_h. The last output index is that starting index plus M + N − 2. MATLAB array positions still begin at one; those positions are storage locations, not necessarily the physical sample indices.

Calculate a Small Example by Hand

Consider x = [1, 2, 3] and h = [1, 1], both starting at sample zero. The first sequence has three samples and the second has two, so we expect four output samples. Writing out the overlap makes the calculation easier to inspect than immediately looking at a graph.

Therefore, y = [1, 3, 5, 3]. Notice that the ending sample has not disappeared just because the input ended. The filter still contributes a final term. Omitting this tail is a common mistake when writing a manual convolution calculator.

If we instead use h = [0.5, 0.5], every output value is halved. We obtain [0.5, 1.5, 2.5, 1.5]. Within the fully overlapping portion, this filter averages two neighboring input samples. At the boundaries, zero padding affects the result, so the first and last values do not represent averages of two available measurements.

Using MATLAB's Built-In Convolution Function

The basic command is y = conv(x,h). MATLAB calculates the full convolution by default. The official conv documentation also describes the 'same' and 'valid' output options. These options select portions of the convolution; they do not turn it into element-by-element multiplication.

x = [1 2 3];
h = [1 1];
y = conv(x,h);          % Expected values: [1 3 5 3]
ySame = conv(x,h,'same');
yValid = conv(x,h,'valid');

The order of the arguments matters when requesting a cropped result. Full convolution is commutative, but a result sized to the first argument need not retain the same samples after swapping arguments. For this tutorial, we use full convolution throughout so that the indexing and all boundary contributions remain visible.

Create the Convolution Calculator in MATLAB

Save the following function as convolution_calculator.m. The file name must match the function name. Place it in your current MATLAB folder or in a folder on the MATLAB path, then call it from the Command Window using the example shown after the code.

This version accepts vectors directly instead of evaluating expressions typed into an interactive input prompt. It validates the input shape, converts arithmetic to double precision and normalizes row and column vectors to the same internal orientation. Matrices and complex vectors are deliberately outside this example's input contract.

function [y,n,maxError] = convolution_calculator(x,h,nx0,nh0)
% Full linear convolution of two finite real vectors.
if nargin < 3
    nx0 = 0;
end
if nargin < 4
    nh0 = 0;
end
validateattributes(x, {'numeric'}, {'vector','real','finite','nonempty'});
validateattributes(h, {'numeric'}, {'vector','real','finite','nonempty'});
validateattributes(nx0, {'numeric'}, {'scalar','real','finite','integer'});
validateattributes(nh0, {'numeric'}, {'scalar','real','finite','integer'});

x = double(x(:).');
h = double(h(:).');
nx0 = double(nx0);
nh0 = double(nh0);
M = numel(x);
N = numel(h);
y = zeros(1,M+N-1);

for i = 1:M
    for j = 1:N
        y(i+j-1) = y(i+j-1) + x(i)*h(j);
    end
end
reference = conv(x,h);
reference = reference(:).';
if any(~isfinite(y)) || any(~isfinite(reference))
    error('Calculation overflowed; reduce the input magnitudes.');
end
maxError = max(abs(y-reference));
nx = nx0 + (0:M-1);
nh = nh0 + (0:N-1);
n = nx0 + nh0 + (0:M+N-2);

fprintf('Manual result:\n');
disp(y);
fprintf('Maximum absolute difference from conv: %.3g\n',maxError);
figure;
subplot(4,1,1); stem(nx,x,'filled'); grid on;
title('Input x[n]'); xlabel('Sample index'); ylabel('Amplitude');
subplot(4,1,2); stem(nh,h,'filled'); grid on;
title('Impulse response h[n]'); xlabel('Sample index'); ylabel('Amplitude');
subplot(4,1,3); stem(n,reference,'filled'); grid on;
title('Full convolution using conv'); xlabel('Sample index'); ylabel('Amplitude');
subplot(4,1,4); stem(n,y,'filled'); grid on;
title('Full convolution using nested loops'); xlabel('Sample index'); ylabel('Amplitude');
end

Now call the function with the hand-calculated example. The expected output values are [1, 3, 5, 3], and the sample indices are [0, 1, 2, 3]. These expected values follow from the arithmetic above; the retained screenshots later in this article show the original version of the tutorial.

[y,n,err] = convolution_calculator([1 2 3],[1 1],0,0);

How the nested loops work

The outer loop selects one sample from x. The inner loop multiplies that sample by every sample in h. Each product belongs at array position i+j-1. We add the product to the value already stored there because several pairs can contribute to the same output sample.

For example, x(2)h(1) and x(1)h(2) both contribute to y(2). The subtraction of one accounts for MATLAB's one-based storage positions. Without it, the first output would be left empty and the final assignment would exceed the intended length.

Preallocating y avoids repeatedly growing the output array. More importantly, it makes the zero-initialized accumulation explicit. Resetting an output sample inside the wrong loop would erase earlier contributions and produce an answer that might look plausible while being mathematically incorrect.

Why compare numerical values as well as plots?

Plots help us understand shape, delay and amplitude, but two similar plots can conceal small numerical differences or an indexing error. The maximum absolute difference checks the actual arrays. An exact integer example with small values should agree exactly in ordinary double-precision arithmetic.

For long sequences containing fractional values, different accumulation orders can introduce small floating-point differences. Judge those differences against the signal scale and the numerical needs of your application. A fixed tolerance suitable for a unit-amplitude signal may be inappropriate for very large numbers or a result dominated by cancellation.

Understand the Output and Original Demonstration

The original Command Window screenshot below shows inputs and the two computed outputs. Its variable names and prompting behavior belong to the archived demonstration. In the revised function, y is the manual result, n contains its sample indices, and maxError measures its difference from conv.

The four-plot arrangement is retained because it lets us inspect both inputs before comparing the outputs. Read the horizontal labels carefully: a curve drawn at array positions one through four is not automatically a signal beginning at sample one. Our function uses the supplied starting indices to make that distinction visible.

The original video remains useful for seeing the convolution workflow and plot layout. Follow the updated function in this article when reproducing the revised interface and validation behavior.

Try These Examples Before Using Longer Signals

1. Use h = [1]. The output should equal x, with its starting index shifted by the starting index assigned to h.
2. Use x = [0, 0, 0]. Every output must be zero, although the full output still has M + N − 1 entries.
3. Use x = [1, 2, 3] and h = [1, −1]. The expected result is [1, 1, 1, −3], which demonstrates signed contributions and the final boundary term.
4. Swap x and h while retaining full convolution. The numerical sequence should remain the same.
5. Supply one input as a column vector. The calculator should return the same row-oriented result.
6. Set nx0 = −2 and nh0 = 3. The first output index should become 1 without changing the output amplitudes.

Work through at least one example with a nonzero starting index on paper. This separates a correct numerical calculation from a correct physical interpretation, which becomes especially important when comparing measured delays or aligning signals from different sources.

Convolution, Correlation and Multiplication

Element-by-element multiplication, written x.*h for compatible arrays, multiplies corresponding entries without summing shifted overlaps. Matrix multiplication, written x*h, obeys matrix-dimension rules and serves a different purpose. Neither is a replacement for full convolution.

Correlation measures similarity as one signal is shifted relative to another. For real sequences, its relationship to convolution involves reversing one input; for complex sequences, conjugation also matters. Lag conventions must be checked before interpreting the sign of a correlation peak. Our speech correlation tutorial uses this idea for a limited template-matching exercise.

Bilinear image interpolation is another weighted calculation, but it should not be described simply as the same correlation algorithm. It estimates pixel values at new spatial coordinates. Understanding these distinctions prevents choosing a familiar command for the wrong mathematical operation.

Sampling, Filtering and Longer Calculations

Sample indices are not automatically seconds

If the sample interval is T_s, convert an index to time using t = nT_s, with any additional physical time origin handled consistently. Both sequences must use compatible sampling intervals for the usual discrete-time system interpretation. Convolving recordings sampled at different rates without resampling does not create a meaningful common time axis.

When approximating a continuous-time convolution integral from uniformly sampled functions, a factor of T_s is normally required: the integral is approximated by a sum of products multiplied by the sample interval. This factor is not automatically required when h already contains discrete-time filter coefficients. Decide which problem you are solving before multiplying the output by a time step.

Moving-average filters and gain

A length-L moving average uses h = ones(1,L)/L. The coefficients sum to one, so a constant input is preserved where the filter fully overlaps it. Omitting the division by L produces a moving sum, whose steady constant-input output is L times larger.

The startup and ending portions of full convolution include zero padding. If you want a streaming filter, also think about stored history between data blocks. Processing each block independently and discarding its tail loses contributions across block boundaries, even though each individual conv call is mathematically correct.

Linear convolution versus circular convolution

Multiplying two discrete Fourier transforms and applying an inverse transform produces circular convolution at the chosen transform length. To recover the full linear convolution of finite sequences, zero-pad both sequences to a transform length of at least M + N − 1 before multiplying their transforms.

A shorter transform allows the tail to wrap around and add into earlier samples. That is time-domain aliasing, not a plotting issue. FFT-based methods can be useful for long signals, but our direct nested loops are intentionally easier to inspect and require M × N multiply-and-accumulate contributions.

Troubleshooting the Calculator

Check a result with independent mathematical properties

For finite sequences, the sum of the full convolution equals the product of the input sums, apart from floating-point error. In our first example, sum(x) is 6 and sum(h) is 2, so the output should sum to 12. The values [1, 3, 5, 3] satisfy that check. This property can expose a missing boundary sample or a dropped contribution without repeating every multiplication.

It is not a complete correctness proof. Two wrong samples can cancel in the total. Combine the sum check with output length, starting index, a known short example and the numerical reference comparison. Each check catches a different kind of mistake and helps explain the calculation rather than merely producing a pass message.

Another useful property is scaling. Multiplying x by a constant should multiply the full result by the same constant. Adding two input sequences should produce the sum of their separate responses through the same h, provided the sequences are aligned and padded consistently. These are the linearity properties that make convolution useful for modeling linear systems.

If you later extend the function to complex data, remove the real-input restriction deliberately and review plotting choices. Convolution itself does not conjugate either input. Plotting only a real component would hide the imaginary component, so provide separate real and imaginary plots or clearly labeled magnitude and phase views.

Review of This MATLAB Convolution Calculator

This calculator is a useful learning tool because it exposes the accumulation, preserves sample origins and provides a built-in reference calculation. Its limitations are equally clear: it handles finite real vectors, uses direct loops, creates plots on every call and does not maintain streaming state.

For a larger project, separate the numerical function from its plotting interface. Save the input parameters with each result, add application-specific tolerances and check representative boundary cases. A graphical interface can then call the same numerical function without duplicating the calculation inside button handlers.

For batch processing, return numerical results without opening a new figure for every signal pair.

Frequently Asked Questions

Can the calculator convolve continuous functions?

It works with numerical vectors. Continuous functions must first be sampled with a clearly chosen interval, and an approximation to the continuous convolution integral needs the appropriate sample-interval factor. Symbolic integration is a different workflow.

Why is the output longer than the input?

Each input sample contributes a shifted copy of the second sequence. The earliest and latest contributions extend the full result to M + N − 1 samples, even when those boundary values happen to be zero.

Does same convolution remove filter delay?

No. It crops the full result to a selected length. You still need to understand the filter's phase, the crop convention and the time origin before interpreting alignment or delay.

Can I use this function for images?

The function intentionally accepts vectors only. Image filtering requires a two-dimensional operation, a kernel and a boundary policy. MATLAB's conv2 is a suitable starting point for studying that separate calculation.

What should I learn after this example?

Try a moving average, a first-difference filter and a delayed impulse. Predict the output before running the calculator, then use both the values and their sample indices to explain any difference from your prediction.