Impulse response synthesis in the frequency domain

Question

Background - Create context

I am trying to synthesise some impulse responses to simulate an active noise control system. I have already formulated the problem in the frequency domain so now I want to move to the time-domain formulation.

The transfer functions I have used in the frequency domain are those of point/monopole sources given by

$$ p \left( \omega, r \right) = \frac{j \omega \rho}{4 \pi r} e^{-j k r} \tag{1} \label{1} $$

where $r$ is the distance between two points (source and receiver), $\rho$ is the density of the medium of propagation (air for my case), $\omega = 2 \pi f$ is the radial frequency with $f$ the temporal frequency, $k = \frac{\omega}{c} = \frac{2 \pi}{\lambda}$ is the wavenumber with $c$ the speed of propagation and $\lambda$ the wavelength and $j = \sqrt{-1}$ is the imaginary unit.

Since I have the transfer function for each frequency of interest from one source to one receiver I would like now to get an impulse response (IR) from this frequency domain representation.

Implementation

To transition to the time domain, I define the following parameters:

1. Sampling frequency: This can be arbitrary based on the needed bandwidth (and possibly other factors).
2. Length of the IR: This will have an impact on the chosen frequencies of interest as it dictates the frequency resolution.
3. Frequencies of interest: These are defined as the frequencies corresponding to the frequency bins in the range of interest (can be a single frequency, a range or discrete frequencies in the range).

To generate the single-sided spectrum I assign the values calculated with equation \eqref{1} in the appropriate bins (whose frequencies are used in the equation to get the transfer functions). Next, I generate the two-sided spectrum by mirroring the first half (except for the DC component) and conjugating. The result for a signal with frequencies from about $10 \, \textrm{Hz}$ to about $80 \, \textrm{Hz}$ can be seen in the image below. The sampling frequency is $300 \, \textrm{Hz}$ and the length of the (double-sided) spectrum vector is $201$ samples. The distance $r$ corresponds to $0.5 \, m$.

All trends expected from equation \eqref{1} can be seen (linear phase and the amplitude increases with frequency). Now, when I take the inverse Fourier transform to acquire the IR of this spectrum I end up with what can be seen in the next figure.

The second plot, termed “Centred” has moved the second half of the IR to the first in a cyclic manner (implementation is performed in MATLAB and uses the ifftshift function). I believe in the second plot, the time axis has to be changed to go from negative half the duration to positive half the duration but this is not done in the figure.

Question

The first, and possibly the most important question here is, whether the second plot of the last figure above, showing the impulse responses, is the correct one. My intuition says that, since the synthesis is performed entirely in the frequency domain, there is no causality constraint and some kind of “time domain aliasing” takes place here but I am not sure I understand it completely.

Furthermore, an additional question is related to the energy spread in time manifesting as ripples in the impulse response. Is there any way to avoid this and concentrate the energy towards the “main part” of the IR?

Answer 1

You can break down your transfer function into two different parts: $\frac{j \omega \rho}{4 \pi r}$ is a differentiator and $e^{-j k r}$ is a pure time delay.

If we were sloppy, we could write the inverse continuous Fourier Transform simply as

$$p(t,r) = \frac{ \rho}{4 \pi r} \delta^{'}(t-r/c) \tag{1} $$

where $\delta^{'}(t)$ is the first derivative of the dirac delta. The sloppiness here comes from the fact that the dirac delta is not a function but a distribution and so is it's derivative.

Sampling this is quite complicated. I would break it into two parts, the differentiator and the delay.

There is a lot of literature on how to implement approximations to a differentiator as discrete filters. For a quick intro you can look here

The delay is easy enough if you can round it to the nearest full sample. If you have to model the delay exactly, then you need to implement a fractional delay, which is also quite complicated. A very thorough treatment of this topic can be found at https://ieeexplore.ieee.org/document/482137 (paper is behind a pay wall but the software is open source).

Is there any way to avoid this and concentrate the energy towards the “main part” of the IR?

You have two main problems here. A differentiator is NOT bandlimited and hence cannot be sampled without aliasing. You need to define what type and amount of aliasing you can tolerate. In your example you have just zero'd some frequency bins. That's equivalent to filtering with a brick wall filter which has a sinc function as an impulse response. That's the main sound source of the time domain ringing and the non-causality.

A second source of ringing would be the delay, if it's not an integer number of samples. Fractional delays are also non-causal.