Tone matching, problem with computation of impulse response

[yvdh]

Hello,

I'm working on a tone-matching impulse response computation tool for acoustic music instruments. The idea is to match the signal generated by the instrument pickup system (typically some kind of piezo-electric device with some less-desirable tonal characteristiscs) with that of a instrument acquired with a good microphone. It's based on the power spectrum of both those signals recorded simultaneously while playing the instrument (we cannot use chirps, white or pink noise here as we cannot 'excite' the instrument using a synthetic signal). By comparing both power spectra, and possibly their noise spectra, a frequency response can be determined. The phases can be reconstructed using either a Hilbert transform approach (not much luck there) or by using a minimum-phase transform approach that removes the out of the unit circle zero's (works in scilab). This is a form of inverse filtering, with all it's numerical instabilities and those are addressed by using a 'Wiener correction', smoothing the spectra and clipping possible gain of the frequency respone curve. I got this working very nicely in SciLab (a free Matlab alternative).
This environment is not great for sharing this tool to a wider audience, hence I am porting this to C# - WPF.

However, I seem to be having issues when reconstructing the impulse response from the computed frequency response in NWaves. Specifically: IR = ifft(exp(fft(Causal(ifft(log(frequency magnitude response))), with causal() consisting of doubling t > 0 values and zeroing t < 0 values.

In NWaves it seems taking ifft of the log spectrum is not working as expected. Looking through the source code it's not very clear to me what the exact format of the input data should be: presence of positive & negative frequencies, order of data ('FFT buffer format' or not, ...).
Any help greatly appreciated, been trying a lot of stuff for 2-3 days without success.
Yves

[ar1st0crat]

Hi, IIUC, this part ifft(log(frequency magnitude response)) is very similar to real cepstrum implemented here .

NWaves provides methods Fft.Direct() / Fft.Inverse() for more "low-level" FFT manipulations:

var fft = new Fft(1024);

re = { ... };  // array of real parts (likely, signal samples), size 1024
im = { ... };  // array of imaginary parts (likely, zeros), size 1024

fft.Direct(re, im);

for (var i = 0; i < re.Length; i++)
{
    var mag = Math.Sqrt(re[i] * re[i] + im[i] * im[i]);

    re[i] = (float)Math.Log(mag + float.Epsilon, logBase);
    im[i] = 0.0f; 
    // or maybe:
    // im[i] = Math.Atan2(im[i], re[i]);
}

fft.Inverse(re, im);

The structure of FFT/IFFT input/output should be identical to default fft()/ifft() in SciLab (symmetric)

If you need double precision, replace float -> double and Fft -> Fft64

[yvdh]

Hi,
thnx for the quick answer.
I knew about the Cepstral transform, but that starts from a time domain signal so it has an extra fft I don't need. I start from a magnitude or power frequency response curve with only the positive frequencies.

I already tried to use the Fft class with all zero phases and by adding the negative frequencies to the arrays in classic fft buffer format (frequency indices {0,1,2,3,..., N/2+1,-(N/2), ...,-1}) but I got very odd results. Mmm, maybe I did something wrong, I will check again.

From a theoretical standpoint it is unclear to me if I should zero the negative frequencies and double the positive one, as well as making the Iff(log(frequency response)) causal by removing the negative time part of that signal. In theory I would think that negative frequencies = negative time so both operations amount to the same thing? But if so the fft(causal in time domaiin (ifft()) operation amounts to just making the negative frequency response frequencies zero ... All a bit confusing, it's been 30 years since I studied DSP (back when the animals were still talking) so I am a bit rusty.

I also knew about the log problems if magnitude is zero (hence the float.epsilon).

I will let you know how it works out!
Yves