Kalibrierung der Ultraschalllautsprecher und Aussendung kalibrierter Signale

Ich versuche, einen Ultraschalllautsprecher mit dem Ziel zu kalibrieren, vorhersagbare Signale zu senden, aber leider habe ich immer wieder Probleme, wahrscheinlich aufgrund meines Mangels an DSP-Fu.

Ein kleiner Hintergrund

Ich möchte in der Lage sein, Sounds so nah wie möglich an einer kalibrierten Aufnahme wiederzugeben, die ich habe. Soweit ich die Theorie verstehe, muss ich die Lautsprecherübertragungsfunktion finden und die Signale, die ich damit aussenden möchte, entfalten. So etwas (im Frequenzbereich):

X -> H -> XH

Wo Xist das ausgesendete Signal Hist die Lautsprecherübertragungsfunktion und XHist XZeiten H. Eine Division ( ./) sollte mir jetzt geben H.

Um nun ein kalibriertes Signal zu senden, sollte es geteilt werden durch H:

X/H -> H -> X

Was wurde getan?

• Platzierter Lautsprecher und ein kalibriertes Mikrofon in einem Abstand von 1 m auf Stativen.
• Aufgezeichnete 30+ lineare Sweeps von 150 kHz bis 20 kHz, 20 ms lang und aufgezeichnet bei 500 KS / s.
• Ausgerichtete und gemittelte Signale mit dem Matlab / Octave-Skript unten, unter dem Skript kann das resultierende Signal angezeigt werden.

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 transformiert Xund XHund die oben genannten Berechnungen durchgeführt, sieht das Ergebnis plausibel aus. Unten ist eine normalisierte Darstellung von H(lila) und X/H(grün).

  

Der Plot wurde auf die relevanten Frequenzen gekürzt.

Bitte lassen Sie mich wissen, wenn ich das falsch mache.

Meine Frage

Nach der Berechnung X/Hmuss ich es zurück in den Zeitbereich transformieren, nahm ich dies ein einfacher wäre , ifft(X./H)und wavwrite, aber alle meine Versuche haben bisher eine plausible Antwort bekommen gescheitert. Ein Frequenzvektor Hf, Hund Xkann gefunden werden hier in MAT7-Binärformat.

Vielleicht bin ich nur müde und es gibt hier eine einfache Lösung, aber im Moment kann ich es nicht sehen. Jede Hilfe / Beratung wird sehr geschätzt.

Ich habe die Antwort gefunden, nachdem ich mir die Referenzen angesehen habe, die Jim Clay in den Kommentaren erwähnt hat. Vielen Dank, Jim.

Ich habe den Fehler gemacht, nur die Größe zu berücksichtigen, die zu einem nullphasigen Signal führt und nicht sinnvoll für die Emission verwendet werden kann, zumindest nicht in diesem Setup.

Der Code, den ich letztendlich verwendet habe, ist unten zu sehen.

Das Skript hält sich an die Namenskonvention, Zeitdomänensignale in Kleinbuchstaben und Frequenzdomänensignale in Großbuchstaben zu halten.

% 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

Die Spektrogramme von x conv hund x deconv hsind unten zu sehen:

Diese erscheinen mir plausibel, obwohl das entfaltete Signal etwas Rauschen enthält.

Der nächste Test wird sein, um zu sehen, ob das Aussenden x_deconv_yetwas Ähnliches ergibt, xabgesehen von den Frequenzen, die der Lautsprecher nicht aussenden kann.

Update mit Testergebnissen

Wir haben die oben beschriebenen Messungen mit einem logarithmischen Down-Sweep überarbeitet. Diese Ergebnisse scheinen darauf hinzudeuten, dass die Methode funktioniert.

Der Verifikationstest bestand aus dem Aussenden X / Hund der Erwartung, Xzurück zu kommen , dh gleiche Energie in allen Frequenzen. Da die schlechteste Ausgangsfrequenz etwa 20 dB schwächer als die beste ist, wird erwartet, dass der höchste Ausgangspegel so viel niedriger ist.

Das Signal, das ausgesendet wurde:

Die Zeitreihen und das Spektrogramm des aufgezeichneten Signals sehen folgendermaßen aus: