Kalibracja głośników ultradźwiękowych i wysyłanie skalibrowanych sygnałów

Próbuję skalibrować głośnik ultradźwiękowy w celu emitowania przewidywalnych sygnałów, ale niestety ciągle mam problemy, prawdopodobnie z powodu braku DSP-fu.

Trochę tła

Chcę móc odtwarzać dźwięki jak najbliżej skalibrowanego nagrania, które mam. O ile rozumiem teorię, muszę znaleźć funkcję przenoszenia głośników i dekonwolować sygnały, które chcę z nią emitować. Coś takiego (w dziedzinie częstotliwości):

X -> H -> XH

Gdzie Xjest emitowany sygnał Hjest funkcja przenoszenia głośników i XHto Xczasy H. Podział ( ./) powinien mi teraz dać H.

Teraz, aby wyemitować skalibrowany sygnał, należy go podzielić przez H:

X/H -> H -> X

Co zostało zrobione

• Umieszczony głośnik i skalibrowany mikrofon w odległości 1 m od statywów.
• Nagrano ponad 30 liniowych przebiegów 150 kHz-20 kHz o długości 20 ms i zarejestrowano przy 500 KS / s.
• Wyrównane i uśrednione sygnały za pomocą skryptu Matlab / Octave poniżej, pod skryptem można zobaczyć wynikowy sygnał.

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 przekształcił się Xi XHwykonał wspomniane powyżej obliczenia, wynik wydaje się wiarygodny. Poniżej znajduje się znormalizowany wykres H(fioletowy) i X/H(zielony).

  

Fabuła została obcięta do odpowiednich częstotliwości.

Daj mi znać, jeśli podejdę do tego w niewłaściwy sposób.

Moje pytanie

Po obliczeniu X/Hmuszę go przekształcić z powrotem w dziedzinę czasu, przyjąłem, że będzie to proste ifft(X./H)i wavwrite, ale wszystkie moje dotychczasowe próby nie dały żadnej wiarygodnej odpowiedzi. A wektor częstotliwości Hf, Hi Xmożna znaleźć tutaj w formacie mat7-binarny.

Może jestem po prostu zmęczony i istnieje tutaj proste rozwiązanie, ale w tej chwili go nie widzę. Każda pomoc / rada jest bardzo mile widziana.

Znalazłem odpowiedź po zapoznaniu się z referencjami, o których wspominał Jim Clay w komentarzach, dziękuję Jim.

Popełniłem błąd, biorąc pod uwagę tylko wielkość, która skutkuje zerowym sygnałem fazowym i nie może być sensownie wykorzystana do emisji, przynajmniej nie w tym ustawieniu.

Kod, którego ostatecznie użyłem, można zobaczyć poniżej.

Skrypt jest zgodny z konwencją nazewnictwa polegającą na przechowywaniu sygnałów w dziedzinie czasu małymi literami i sygnałów w dziedzinie częstotliwości dużymi literami.

% 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

Spektrogramy x conv hi x deconv hmożna zobaczyć poniżej:

Wydaje mi się to prawdopodobne, chociaż w rozłożonym sygnale występuje pewien szum.

Następnym testem będzie sprawdzenie, czy emitowanie x_deconv_ydaje coś podobnego z xwyjątkiem częstotliwości, których głośnik nie może emitować.

Zaktualizuj z wynikami testu

Ponownie zmierzyliśmy opisane powyżej pomiary, używając logarytmicznego przemiatania w dół. Te wyniki wydają się sugerować, że metoda działa.

Test weryfikacyjny polegał na emisji X / Hi oczekiwaniu na Xpowrót, tj. Równej energii na wszystkich częstotliwościach. Ponieważ najgorsza częstotliwość wyjściowa jest o około 20 dB słabsza niż najlepsza, najwyższy poziom wyjściowy powinien być znacznie niższy.

Sygnał, który został wysłany:

Szeregi czasowe i spektrogram zarejestrowanego sygnału wyglądają następująco: