What's the correct syntax to set-up discrete convolution ?

I am trying to generate a 1-D discrete convolution data-set but I keep getting NaN values. I am wondering what's going wrong. I want to convolve the function "differentiated_green" with the function "f" with respect to the variable "x".

These are the common parameters needed for complete definition of both the functions:

M = 0.5; alpha = log(2)/2; epsilon = 0.5; f = 1/30; omega = 2*pi*f; c_inf = 340; k = omega/c_inf; t = 270; nodes = 401;
y = 0;
x_start = -200; x_end = 200;
length = x_end - x_start;
x = linspace(x_start,x_end,nodes);

This is the definition of the function "f":

f = epsilon * exp( -alpha * x.^2);

This is the definition of the function "differentiated_green":

A = omega/( 4 * (c_inf^2)*sqrt(1-M^2));
B = besselh(0, (k*x)/(1-M^2));
C = exp( -1i*((M*k*x)/(1-M^2)) - (1i*omega*t));
differentiated_Green = (A * B) .* C;

The function "differentiated_green" is a complex value due to the presence of Hankel function in "B" and the exponential term "C". Since I want to convolve this function's real values with the function "f", these are the commands I am using:

p = conv(real(differentiated_Green),f,'same');

These commands are giving me NaN values and I am trying to understand what's going wrong. Can someone please point out to me where I should troubleshoot ?

Answer 1

[Paul]

The variable B contains one or more NaN values. Use isnan to find which ones and then review the the corresponding inputs to besselh and see if those are NaN or, if not, figure out why besselh returns NaN for those inputs.

[Walter Roberson]

find(isnan(B)) shows 201 -- corresponding to x == 0

[ishan]

Apologies for the late reply but I am still struggling with this. There is a NaN value at the centre of the domain and it's getting carried everywhere. Can you suggest any toy problems which involve convolution with a Hankel function ? Would there be any analytical solution for comparison ? I want to further my understanding of conv and conv2 functions so that i can use them properly.

[ishan]

To add, my verification of 1D and 2D convolution of two Gaussians is working perfectly well because I can compare the result with analytical expressions.

[Walter Roberson]

M = 0.5; alpha = log(2)/2; epsilon = 0.5; f = 1/30; omega = 2*pi*f; c_inf = 340; k = omega/c_inf; t = 270; nodes = 401;
y = 0;
x_start = -200; x_end = 200;
length = x_end - x_start;
syms x
f = epsilon * exp( -alpha * x.^2);
A = omega/( 4 * (c_inf^2)*sqrt(1-M^2));
B = besselh(0, (k*x)/(1-M^2));
C = exp( -1i*((M*k*x)/(1-M^2)) - (1i*omega*t));
differentiated_Green = (A * B) .* C

differentiated_Green = 

limit(real(differentiated_Green), x, 0, 'left')

ans = 

limit(real(differentiated_Green), x, 0, 'right')

ans = 

The limits are negative of each other.

So x == 0 has a non-removable discontinuity. Therefore it makes no sense to convolve the function over an interval that includes the discontinuity.

[Paul]

@ishan

Is conv being used truly for a discrete convolution sum? Or is it being used to approximate a continuous convolution integral?

[ishan]

@Walter Roberson Thanks for laying this out. I should have realized this right from the start. @Paul When you say "continuous" are you asking if there is an analytical solution to the convolution ? I don't know to be honest. Sorry, but I did not understand when you asked whether conv is truly being used for discrete convolution.

[Paul]

By "continuous" I mean that the domain of the functions to be convolved is the entire real line. In this case, the convolution of f(x) and g(x) is given by the continuous convolution integral and the domain of the result is the entire real line.

If the domain of the functions f[n] and g[n] to be convolved are the integers, then their convolution is given by the discrete convolution sum and the domain of the result is the integers.

In many (all?) cases, the convolution sum applied to uniformly spaced samples of f(x) and g(x) can be used to approximate uniformly spaced samples of the convolution integral of f(x) and g(x).

As I understand, the domain of the Hankel function is the entire real line (actually the complex plane), with the possible exception of the origin. If we assign a value to the Hankel function at the origin, then we might (I'm not sure) be able to define a convolution, but the first thing is to determine if you want the convolution integral or the convolution sum.