Ultrasonic speaker calibration and emitting calibrated signals

Question

I'm trying to calibrate an ultrasonic speaker with the aim of emitting predictable signals, but alas I keep running into trouble, probably due to my lack of DSP-fu.

A little background

I want to be able to playback sounds as close as possible to a calibrated recording I have. As far as I understand the theory, I need to find the speakers transfer function and deconvolve the signals I want to emit with it. Something like this (in the frequency domain):

X -> H -> XH

Where X is the emitted signal H is the speakers transfer function and XH is X times H. A division (./) should now give me H.

Now, in order to emit a calibrated signal, it should be divided by H:

X/H -> H -> X

What's been done

• Placed speaker and a calibrated microphone 1 m apart on tripods.
• Recorded 30+ linear sweeps 150KHz-20KHz, 20ms long, and recorded @ 500 KS/s.
• Aligned and averaged signals with the Matlab/Octave script below, under the script the resulting signal can be seen.

files = dir('Mandag*');

rng = [1.5e6, 1.52e6];

[X, fs] = wavread(files(1).name, rng);
X = X(:,1);

for i=2:length(files)
    [Y, fs] = wavread(files(i).name, rng);
    sig = Y(:,1);
    [x, off] = max(xcorr(X', sig'));
    off = length(X) - off;
    if(off < 0)
        sig = [zeros(1, -off), sig(1:end+off)'];
    elseif (off > 0)
        sig = [sig(off:end)', zeros(1, off-1)];
    end
    X = X + sig';
end
X = X/length(files);

• Fourier transformed X and XH and did the calculations mentioned above, the result looks plausible. Below is a normalized plot of H (purple) and X/H (green).

  

The plot has been truncated to the relevant frequencies.

Please let me know if I'm going about this the wrong way.

My question

After calculating X/H I need to transform it back to the time domain, I assumed this would be a simple ifft(X./H) and wavwrite, but all my attempts so far have failed to get any plausible answer. A frequency vector Hf, H, and X can be found here in mat7-binary format.

Maybe I'm just tired and there is a simple solution here, but at the moment I can't see it. Any help/advice is much appreciated.

Answer 1

Found the answer after looking at the references that Jim Clay mentioned in the comments, thank you Jim.

I made the mistake of only considering the magnitude which results in zero-phased signal and can not sensibly be used for emission, at least not in this setup.

The code I finally ended up using can be seen below.

The script adheres to the naming convention of keeping time domain signals in lower-case and frequency domain signals in upper-case.

% Align and sum all files called Mandag*
files = dir('Mandag*');

% Where in the recordings the signal is
rng = [1.5e6, 1.52e6];

% Initialize the xh vector
[xh, fs] = wavread(files(1).name, rng);
xh = xh(:,1);

for i=2:length(files)
    y = wavread(files(i).name, rng);
    y = y(:,1);
    % Determine offset between xh and y
    [~, off] = max(xcorr(xh', y'));
    off = length(xh) - off;
    % Shift signal appropriately
    if(off < 0)
        y = [zeros(1, -off), y(1:end+off)'];
    elseif (off > 0)
        y = [y(off:end)', zeros(1, off-1)];
    end
    xh = xh + y';
end

% Average
xh = xh/length(files);

% Location of the 20ms signal
xh = xh(2306:12306-1);

% Normalize
xh = xh / max(xh);

% Apply a moving average filter on xh to reduce noise. Window size of 4 was
% experimentally determined to give the best results
n = 4;
B = zeros(n, 1);
for i=1:n
  B(i) = 1/n;
end
xh = filter(B, 1, xh);
xh = xh / max(xh);

x = wavread('sweep.wav');
x = x(1:2:end);            % Sweep generated @ 1MHz, decimate
                           % to have same length as xh

% Transform x into frequency domain and determine H
X = fft(x);
H = fft(xh) ./ X;

% Vector indices to choose only frequencies of interest
starti =  20e3 / 50;
endi   = 100e3 / 50;
rng    = starti:endi;
irng   = (length(x) - endi) : (length(x) - starti);

% Zero out unwanted frequencies
X = [zeros(1,      starti - 1   ), X( rng)', zeros(1, length(X)/2 - endi) ...
     zeros(1, length(X)/2 - endi), X(irng)', zeros(1,      starti - 1   )]';

% Deconvolve x with h
X_deconv_H = X ./ H;

% Transform X/H to time domain and normalise
x_deconv_h = real(ifft(X_deconv_H));
x_deconv_h = x_deconv_h / max(x_deconv_h);

% Save the deconvolved sweep
wavwrite(x_deconv_h, fs, 'deconvolved_sweep.wav');

% Generate  spectrograms of xh and x_deconv_h
winsize = 512;
overlap = round(.99 * winsize);
figure(1)
specgram(xh, winsize, fs, hann(winsize), overlap)
colorbar
figure(2)
specgram(x_deconv_h, winsize, fs, hann(winsize), overlap)
colorbar

The spectrograms of x conv h and x deconv h can be seen below:

These seem plausible to me although there is some noise in the deconvolved signal.

Next test will be to see if emitting x_deconv_y gives something resembling x barring those frequencies that the speaker cannot emit.

Update with test results

We redid the measurements described above using a logarithmic down-sweep. These results seem to suggest that the method works.

The verification test consisted of emitting X / H and expecting to get X back, i.e. equal energy in all frequencies. As the worst output frequency is about 20dB weaker than the best, the highest output level is expected to be that much lower.

The signal that was emitted:

The time series and spectrogram of the recorded signal look like this: