Volterra Kernel Convolution Method

Question

Here I understand that the first term is a simple convolution of the input signal with first volterra kernel that I have acquired through farina sweep method, my question is if the second term is the convolution of second volterra kernel with our input signal raised to the power of two? Would the code using matlab, of the first and second terms be correctly written as `

outputsignal1=conv(voltkernel1,inputsignal.^1); %basic convolution with system impulse
outputsignal2=conv(voltkernel2,inputsignal.^2); %convolution of second kernel with squared signal
summedsignal=outputsignal1+outputsignal2;  %dynamic nonlinear output signal with memory

Answer 1

Yes the MATLAB code is correct. Be careful though, the bandwidth of the signal squared is twice that of the signal itself, which may lead to aliasing if the sampling frequency is too low compared to the signal bandwidth. This can be remedied by properly resampling the signal to a higher sampling frequency with a lowpass filter, then squaring, low-pass filtering, and decimating back to the original rate.

Another pitfall if you are doing everything by yourself is that the impulse-response-like partitions of the deconvolved sweep response are not directly the Volterra kernels. The rest of this answer will discuss how to get usable kernels from the partitions corresponding to the harmonics.

First, you may need to phase shift each partition with a constant phase shift, depending on what kind of sweep you used. You can test if there is any phase shift needed (and test if you do the phase shift correctly), by using as a test system a memoryless non-linearity that is the sum of a number of early Chebyshev polynomials of the first kind, as this will generate for a unit amplitude cosine input a number of harmonics in the same phase (explanation further down in this answer). If I understood the following paper correctly, for certain kinds of sweeps there is no need for phase shifts:

Antonín Novák, Laurent Simon, František Kadlec, Pierrick Lotton. Nonlinear System Identification Using Exponential Swept-Sine Signal. IEEE Transactions on Instrumentation and Measurement Volume: 59, Issue: 8, Aug. 2010.

Then each one-dimensional Volterra kernel can be calculated as the weighted sum of the (phase-shifted) partitions, with the weights obtained from solving a simple linear system. This is covered in Angelo Farina's papers, but at least some of them complicate the derivation by using sines rather than cosines, which gets into the earlier phase shift territory. Cosines are simpler because they have real frequency-domain representations implying no phase shift. Here is an Octave script (using cosines) that solves the weights given the number of harmonics:

N = 3; #Number of harmonics including the fundamental (you can touch this)
M = 2*N+1; #Matrix height, must fit positive, zero and negative frequencies
x = zeros(M, N+1);
# From cosines to complex exponentials. Negative indexing by shifted columns
x(1, 1) = 1; #DC bias
x(3, 2) = 0.5; #This is negative fundamental
x(1, 2) = 0.5; #This is positive fundamental
for k = 3:N+1
  x(:, k) = conv(x(:, k-1), x(:, 2))(1:M); #Time domain exponent k
endfor
for k = 1:N+1
  x(:, k) = shift(x(:, k), 1-k); #Shift back columns
endfor
x = x(1:N+1, :); #Go back to cosines
x(2:N+1,:) .*= 2 #
inv(x) #Solve system

The output is:

x =
   1.00000   0.00000   0.50000   0.00000
   0.00000   1.00000   0.00000   0.75000
   0.00000   0.00000   0.50000   0.00000
   0.00000   0.00000   0.00000   0.25000

ans =
   1  -0  -1  -0
   0   1  -0  -3
   0   0   2  -0
   0   0   0   4

The answer ans is the inverse of matrix x. The answer means that if you have three harmonics including the fundamental, then you go from the partitions $h_n$ corresponding to harmonics $n \in \{1,\, 2,\, 3\}$ to kernels $k_m$ corresponding to exponents $m \in \{1,\, 2,\, 3\}$ by:

$$\left\{\begin{eqnarray}k_0 &=& h_0-h_2\\ k_1 &=& h_1 -3h_3\\ k_2 &=& 2 h_2\\ k_3 &=& 4 h_3\end{eqnarray}\right.$$

This is the solution of:

$$\begin{align}&k_0 + k_1 \cos(t) + k_2 \cos^2(t) + k_3 \cos^3(t) = h_0 + h_1 \cos(t) + h_2 \cos(2t) + h_3 \cos (3t)\\ \Rightarrow & k_0 + \frac{k_2}{2} - h_0 + \left(k_1 + \frac{3k_3}{4} -h_1\right)\cos(t) + \left(\frac{k_2}{2} - h_2\right)\cos(2t) + \left(\frac{k_3}{4} - h_3\right)\cos(3t) = 0,\end{align}$$

where $h_0$ is a zero frequency bias that is not determined by the sweep method and is normally assumed zero, $h_0 = 0,$ and $k_0$ is a term that eliminates the bias. I'm not sure how to use it, probably just add to each output sample the sum of samples of $k_0.$

The solution matrices for smaller N are contained in the solution matrices for larger N. The matrices follow the On-line Encyclopedia of Integer Sequences (OEIS) sequence A053120. OEIS gives simpler ways to calculate the weights, and also hints to a link to Chebyshev polynomials of the first kind. They have coefficients that are identical to the weights. The link is actually such that if you use as the memoryless non-linear elements not powers but Chebyshev polynomials of the first kind, $T_n(x)$, then the kernels to use are directly the (phase-shifted) partials corresponding to the harmonics. This comes from the property $T_n\big(\cos(t)\big)=\cos(nt).$ This is also covered in:

Antonin Novak, Laurent Simon, Pierrick Lotton, Joel Gilbert. Chebyshev Model and Synchronized Swept Sine Method in Nonlinear Audio Effect Modeling. Proc. of the 13th Int. Conference on Digital Audio Effects (DAFx-10), Graz, Austria, September 6-10, 2010.

The N = 7 solution from the script agrees with:

Marc R´ebillat, Romain Hennequin, Etienne Corteel, Brian Katz. Prediction of harmonic distortion generated by electro-dynamic loudspeakers using cascade of Hammerstein models. 128th Convention of the Audio Engineering Society, May 2010, London, United Kingdom. pp.7993, 2010. <hal-00619357>

Answer 2

I think actually your code is not correctly implementing the equation you gave. What you are actually imeplementing in matlab is this:

$$ y(n) = \sum_{m_1=0}^{M-1} h_1(m_1)x(n-m_1) + \sum_{m_1=0}^{M-1}h_2(m_1)x(n-m_1)x(n-m_1) $$

If you are getting your higher order kernels from a sine sweep, this is what is available to you, as your higher order "impulse responses" are still one-dimensional.

The model you are using here is not the "full" volterra model, but a simplified model that I have heard referred to as a Hammerstein model in the past.

For the "full" Volterra model, higher order impulse repsonses are N-dimensional, where N is the order.

Answer 3

As some users have told, it's important to stress that the non-linear convolution result gives you (once windowed) the high frequency responses, not the Generalized Hammerstein/Volterra Diagonal ones. You have to convert them.

Another important thing often not explicitly mentioned in papers and that gave me a lot of headaches:

• the maximum f2 is f2 = fs/(2*N) you need enough Nyquist frequency (fs/2) to contain the harmonics of f2

• the minimum frequency you can simulate is f1*N

The brilliant post of Olli Niemitalo includes an easy workaround about the high fs sampling needed that makes me think the classic "why didn't I figure out this until now?"

Matlab resources:

1. Hammerstein Toolbox: Novak (2010) and Rebillat (2011) implementations. No fade-in and fade-out https://www.mathworks.com/matlabcentral/fileexchange/30897-hammerstein-toolbox

2. Novak's site implementation (2010) http://ant-novak.com/Swept_Sine_Simulation.php

3. Schmitz's HKISS/OHKISS toolbox (2016) (available in GitHub) Sadly buggy, with an alternative sine sweep definition and features to actually simulate the resulting GHammerstein model