How to correct phase of Diagonal Volterra kernels obtained with exponential sweep sine?

Question

I'm trying to model a non-linear system using non-linear convolution with Novak's (2010) synchronized exponential sine weep (SESS) that models them with a Generalized Hammerstein (Volterra diagonal).

I get pretty clear kernels in terms of SNR but clearly phase-distorted. It's a common problem in real estimation.

I can theoretically correct them with a time compacting inverse Kirkeby filter (Farina, 2007).

[EDITED AFTER FIRST ANSWER] I use Novak's himself code from http://ant-novak.com/Swept_Sine_Simulation.php

Changes: Doubled both number of kernels and sampling frequency and raised a lot SESS duration for better kernels separation.

fs = 192000;         % sampling frequency
f1 = 10;            % start frequency
f2 = 12000;         % final frequency
T_ = 50;            % approximative time duration
N = 8;              % order of nonlinearities

L = round(f1/log(f2/f1)*T_)/f1;
T = L*log(f2/f1); % this definitive sweep duration in theory makes Novak's SESS more robust to phase issues (synchronized)
t = (0:ceil(fs*T)-1)./fs;
x = sin(2*pi*f1*L*(exp(t/L)-1)); % sweep

Novak uses an elegant frequency domain definition alternative of the regular inverse sweep with amplitude correction. I've tested both with similar results, although Novak's frequency method gives better SNR.

The other change is that I model a "real" system, well, only a guitar amplifier VST plugin. So I use the plugin's output to the SESS for the deconvolution stage.

I'm not 100% sure, but the kernels look extremely "right-sided" like in Farina's 2007 paper. The frequency magnitude of the kernels is coherent with a high gain guitar amplifier VST (even harmonics stronger). The kernels phase is a mess, triangle-saw like. When I test the model with a dry guitar sound as input the output barely improves a regular (first kernel) equalization matching.

Both in the SESS (plugin input/output) and simulation (dry guitar signal) stages I apply a steep bandpass in the range 100hz-10Khz. The plugin is configured to do x2 oversampling (384Khz) to avoid aliasing artifacts.

Answer 1

A potential pitfall in non-linear modeling is using different signal levels. The sweep test needs to be carried out at a sufficiently high signal level to cover the maximum possible real-world input levels, but not much higher as this will require more kernels for the same effective quality. Usually real-world signals to the system are much of the time at lower level than their peak level, with the 1st order linear kernel dominating.

To me your kernels look just fine. I don't know how common it is to see that a higher-order kernel (3rd order in your plots) is more prominent than the 1st order kernel. The saw-like phase comes from the large initial delay, which gives a steep linear dependency of phase on frequency. Wrapped phase of $\pi$ is the same as a wrapped phase of $-\pi$, and you are seeing wrapped phase, which gives the sawtooth pattern. If you remove the initial delay in the inspection of the results, you should be able to unwrap the phase more easily for inspection.

Most systems you would model by the technique will have right-sided kernels. The systems are causal and cannot store energy for very long after they have received an impulse as input. That would require a delay line. Instead, much of the energy passes through almost instantly as the initial spike, and the residual energy is released as the tail.

I suggest you test with a system you have complete control over, for example a memoryless non-linearity that is the sum of the $N$ first Chebyshev polynomials of the first kind, one for each of the $N$ harmonics you want to generate, including the 0th and 1st harmonics in the count. Then any phase shifts due to errors in the algorithms/software should be clearly visible. If the algorithms work correctly, each harmonic impulse response should be a delayed sinc, and the the coefficients of the memoryless nonlinearity polynomial should be resolved as the amplitudes of the scaled impulse-like kernels.

What I have done in my own implementation is to generate the sweep and its fade in and fade out in the frequency domain by programming the timing into the group delay, and amplitude into the magnitudes, rather than using a direct time domain equation for the sweep. Either way should work and will have its own problems. The frequency-domain approach gives amplitude oscillations near the fades. In the frequency domain it is easy to correct the non-whiteness of the sweep in the inverse filter. (I can't publish my implementation in detail at the moment.)

If there are phase problems, I'd try $\cos$ instead of $\sin$ in the time-domain sweep formula, and fade the sweep in and out. The fades will take care of eliminating transients at the ends. The initial and final phases of the sweeped sinusoid could be made such that the value of the sweep is zero both in the beginning and in the end, but with the fades in and out that is no longer needed for transient reduction, because the fades gently bring the value to zero in any case.

I don't see a problem with the sweep going all the way to the Nyquist frequency, if the non-linearity is analog and the analog-to-digital converter (ADC) has an anti-alias low-pass filter at the Nyquist frequency. This way the harmonic impulse responses and consequently the Hammerstein kernels have only non-aliased frequencies up to Nyquist frequency. It is sufficient information to model the non-linear system in the sub- Nyquist frequency band if the system is compatible with the Hammerstein model. Only the implementation has to be careful not to introduce aliasing. In any case, you'd typically want the linear kernel to have frequency response information up to or at least near the Nyquist frequency, which dictates the sweep end frequency.

Answer 2

Clearly I'm missing something. One of the biggest sources of artifacts is the last sample of the SESS not being zero. A fade-out needs to be applied.

In the most important Novak's paper (2015) the definition of the SESS is changed removing the "-1" element of the equation. Sadly, it's a paid paper I don't have. I only know it from a posterior paper (2016), in which Novak states:

L = T/log(f2/f1);
x = sin(2*pi*f1*L*(exp(t/L)));

Quoting Novak:

"Note that the definition of the exponential swept-sine does not contain the " − 1" term contrary to the usual definition [Novak 2010]."

"It has been shown, in [Novak 2015], that in removing the " − 1" term from the original definition, the synchronized swept-sine becomes suitable for estimation of the HHFRs without any need of correction or restriction of the parameters f 1 , f 2 , and T."

This means no sweep synchronization, the kernels are no longer phase-corrected. Freedom to choose the parameters f1, f2 and length. I also guess that freedom of fs as long f2 is below Nyquist frequency and no fade-out. Excellent, but this should have major problems like mixing fundamental and harmonic frequencies.

Probably that's what the latest author about the subject (Schmitz, 2016) uses. Schmitz developed a complete Matlab toolbox (OHKISS) using an ESS until almost the Nyquist frequency while obtaining at least 10 kernels (f2 = 22049, fs = 44100).

Sadly the OHKISS code is buggy and, again, the important paper with the proper math explanation is not available. It also uses advanced GUI features I'm not able to port to Octave.

Answer 3

Many thanks @OlliNiemitalo. I'm learning a lot and hope to resolve my problem soon thanks to you.

I'm aware of the problems of input level dependency and I see you're working on it. Dynamic convolution, like the one making money for the Acustica Audio developers. I'm still far from that level. Just a few doubts and comments:

• Definitely there's a problem with the higher order harmonics. I'm almost sure there's a problem with the fade-out that boosts a lot certain harmonic kernels. I think my simulations doesn't sound good because this overpowered harmonics are masking the effect of the proper or first ones. That may explain why I get too much highs.

• What are your thoughts about the implementations using sweeps with f2 close to the Nyquist frequency? Novak's 2015 seems to use it, others as well. With the classic one I need a very high sampling rate to compute more than 4 kernels with limited bandwidth.

• There seems to be a very desirable property: that the amplitude of the final sample of the ESS be it zero (zero phase depending of the choice of ESS). The first one too, but it seems to give far less artifacts. Schmitz uses a Mixed NonLinear Integer Program solver (SCIP) only for this, optimizing f1, f2 and length given a range.

• Is your frequency domain design similar to the Novak's one when computing the inverse amplitude corrected sweep? At first sight it doesn't seem like to me.

[Added after more researching] Although my problems in the test stage aren't resolved yet, I've discovered the famous implementation of the SESS design in [Novak 2015] in Di Rosario's 2011 paper. This method indeed achieves zero phase (zero amplitude) both in the convenient points for phase corrected kernels and in the end of the SESS. That gives the lowest high end frequency ripple (0.8dB) without using fade-out. Fade-in seems to be mandatory to avoid big low end frequency ripples, no matter what start phase the SESS has.

When I have time (sigh) I'll try to find the mistake in my implementation, which could be perhaps just the weighting of the different kernels. That'd explain the good appearance of the IRs both in magnitude and phase shapes (the latter thanks to Olli explanation) but too high magnitudes in higher harmonics.

I hope the time I'm spending myself and the ones who kindly answer will help someone )