Take the 2-minute tour ×
Signal Processing Stack Exchange is a question and answer site for practitioners of the art and science of signal, image and video processing. It's 100% free, no registration required.

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);

Aligned and averaged signals

  • 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).

    Frequency plot of H and X/H

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.

share|improve this question
1  
This thread - dsp.stackexchange.com/questions/953/… - and this one- dsp.stackexchange.com/questions/2705/… - might be useful to you. –  Jim Clay Jul 20 '12 at 19:23
    
Yes, found my mistake thank you Jim. I was only considering the magnitude of the signals, which results in a zero-phased time signal. Seems like I've got the right result now and I'll add that as an answer. –  Thor Jul 20 '12 at 23:21
add comment

1 Answer

up vote 3 down vote accepted

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:

spectrogram of x conv h spectrogram of x deconv h

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:

Time series of emitted signal

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

Time series of emitted signal Time series of emitted signal

share|improve this answer
    
Any update on this? How did you emitted signal look? –  Busk Aug 1 '12 at 14:44
    
@Busk: Thanks for the interest. I haven't had a chance to test it yet as the equipment is being used elsewhere. I'll post the results when I've done the test. –  Thor Aug 1 '12 at 16:04
    
@Busk: we've now tested it and it seems to work :-). –  Thor Oct 31 '12 at 8:33
    
Looks good! Thanks for the update –  Busk Nov 1 '12 at 5:31
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.