Convolution of function on a non symmetric axis using 'conv'

Hi all,

I'm struggling with a little convolution problem. I have to fit some data with a function that is the convolution of a function R(x) with a certain model M(x). I have noticed that when my data are centered simmetrically around zero all is fine but as soon they are shifted from zero the convolution calculated using 'conv' shifts away of the same amount (see code and example below).

The result does not differs if I modify my M function with a different one and by changing the value of x0 (center of the model function and R(x)) also the convolution with the symmetric x-axis shifts.

Since I'm using symmetric functions I didn't expect this behaviour. Any ideas?

close all; clear all; clc

x = -0.2:1e-4:0.2;  % x-axis
xx = x+0.01;        % shifted x-axis (non symmetric)
s = 1e-2; x0 = 0;   % model parameters

% define a resolution function
R1=exp(-((x-x0)./1e-4).^2); R1=R1/max(R1);
R2=exp(-((xx-x0)./1e-4).^2); R2=R2/max(R2);

% define a model function
G1=exp(-((x-x0)./s).^2);
G2=exp(-((xx-x0)./s).^2);

d1 = diracF(x,x0); % delta
d2 = diracF(xx,x0); % delta

M1 = .7*d1 + 0.3*G1;
M2 = .7*d2 + 0.3*G2;

% convolve model with resolution
A1 = conv(M1,R1,'same');
A2 = conv(M2,R2,'same');

figure; 
subplot(2,1,1); hold on; set(gca,'YScale','log'); ylim([1e-3 2]); xlim([min(x) max(x)])
    plot(x, R1/max(R1), 'g--', 'LineWidth', 1);
    plot(x, M1/max(M1), 'r--', 'LineWidth', 1);
    plot(x,A1/max(A1),'k')

subplot(2,1,2); hold on; set(gca,'YScale','log'); ylim([1e-3 2]); xlim([min(xx) max(xx)])
    plot(xx, R2/max(R2), 'g--', 'LineWidth', 1);
    plot(xx, M2/max(M2), 'r--', 'LineWidth', 1);
    plot(xx,A2/max(A2),'k')
    
legend({'$R(x)$','M(x)','$ R(x)*M(x)$'},'Interpreter','latex','Box','off')
xlabel('x')

function y = diracF(t,t0)
    [~,closestIndex] = min(abs(t-t0));
    y = zeros(1,max(size(t)));
    y(closestIndex) = 1;
end

[Matt J]

Shouldn't R1 be defined as,

R1=exp(-((x-x0)./1e-4).^2); R1=R1/max(R1);

[Matt J]

Since I'm using symmetric functions I didn't expect this behaviour

Why not? It is a well-known property of convolution that if you time-shift one of its operands, the convolution result also shifts. That's why convolutions are used to describe time-invariant systems.

[MMG1]

Thanks. Yes, you're right. And yes, R1 is actually defined on x and not xx, it's just a typo. In my case the point is that both my R and M functions are not on a symmetric range and I don't have symmetric data. So, I guess I'll have to add dummy points to my real data in order to have a symmetric range before doing any fitting.

Answer 1

[Matt J]

This might be what you're after:

% convolve model with resolution

A1 = axconv(M1,R1,x); 
A2 = axconv(M2,R2,xx); 


figure; 
subplot(2,1,1); hold on; set(gca,'YScale','log'); ylim([1e-3 2]); xlim([min(x) max(x)])
    plot(x, R1/max(R1), 'g--', 'LineWidth', 1);
    plot(x, M1/max(M1), 'r--', 'LineWidth', 1);
    plot(x,A1/max(A1),'k')

subplot(2,1,2); hold on; set(gca,'YScale','log'); ylim([1e-3 2]); xlim([min(xx) max(xx)])
    plot(xx, R2/max(R2), 'g--', 'LineWidth', 1);
    plot(xx, M2/max(M2), 'r--', 'LineWidth', 1);
    plot(xx,A2/max(A2),'k')
    
legend({'$R(x)$','M(x)','$ R(x)*M(x)$'},'Interpreter','latex','Box','off')
xlabel('x')

function y = diracF(t,t0)
    [~,closestIndex] = min(abs(t-t0));
    y = zeros(1,max(size(t)));
    y(closestIndex) = 1;
end

function A=axconv(M,R,t)
%Convolution with a specified t axis that includes t=0

 abst=abs(t);
 t0=find(abst==min(abst));

 A=circshift(conv(M,R,'full'),1-t0);
 A(numel(t)+1:end)=[];
   
 
end

[MMG1]

Many thanks!

Answer 2

[Paul]

Hi MMG1

It looks like the code is using the convolution sum to approximate a convolution integral. Is that correct? If so, my comments are:

The Symbolic Math Toolbox can be used to compute a closed form expression for the convolution integral for these examples.

More importantly, M1 is the weighted sum of a dirac and and an exponential. Approximating a convolution integral using conv generally requires scaling the ouput of conv by dx. One exception is when convolving with a dirac, in which case the dirac is replaced with a discrete unit pulse, as in diracF, and no scaling of conv is required. So, in this case I think we need to call conv twice, once with 0.3*G1 and the other with 0.7*d1, with the output of the former scaled by dx before adding them together to get the final result.

dx = 1e-4;
x = -0.2:dx:0.2;  % x-axis

s = 1e-2; x0 = 0;   % model parameters

% define a resolution function
R1 = exp(-((x-x0)./1e-4).^2); R1=R1/max(R1);

% define a model function
G1 = exp(-((x-x0)./s).^2);

d1 = diracF(x,x0); % delta

M1 = .7*d1 + 0.3*G1;

Convolve the two pieces separately. Scale the G1 piece by dx to approximate the convolution integral. Do not scale the d1 piece.

% convolve model with resolution
A1 = conv(M1,R1,'same');
A1g = conv(0.3*G1,R1,'same')*dx;
A1d = conv(0.7*d1,R1,'same');

figure;
ax1 = subplot(211);
hold on;
plot(x, A1g, 'k')

ax2 = subplot(212);
hold on
plot(x, A1d, 'k')

clear R1 G1 M1 A1 A1g A1d
syms x real
s = sym(1e-2); x0 = 0;   % model parameters
% define a resolution function
R1(x) = exp(-((x-x0)./sym(1e-4)).^2); % R1(x) = R1/max(R1)
R1(x) = R1(x)/R1(x0);

G1(x) = exp(-((x-x0)./s).^2);

M1(x) = .7*dirac(x-x0) + 0.3*G1(x);

syms w real
A1(x) = int(M1(w)*R1(x-w),w,-inf,inf)

A1(x) = 

Two separate convolution integrals to compare with the results from conv

A1g(x) = int(0.3*G1(w)*R1(x-w),w,-inf,inf);
A1d(x) = int(0.7*dirac(w-x0)*R1(x-w),w,-inf,inf);

Overplot.

fplot(ax1,A1g(x),[-0.1 0.1]);
fplot(ax2,A1d(x),[-0.1 0.1]);
xlim(ax2,[-0.5 0.5]*1e-3);

The dx spacing might be a little large based on the bottom plot, I guess that depends on what you need.

function y = diracF(t,t0)
    [~,closestIndex] = min(abs(t-t0));
    y = zeros(1,max(size(t)));
    y(closestIndex) = 1;
end