我正在尝试校准超声波扬声器,以发出可预测的信号,但可惜我一直遇到麻烦,这可能是由于缺乏DSP-fu造成的。

背景有点儿

我希望能够播放尽可能接近我已校准的录音的声音。据我了解的理论,我需要找到扬声器的传递函数,并解卷积我想随其发出的信号。像这样(在频域中):

X -> H -> XH


其中X是发射信号H是扬声器传递函数,而XHX乘以H。现在,除法(./)应该给我H

现在,要发出校准信号,应将其除以H

X/H -> H -> X


完成的操作


在三脚架上相距1 m放置了扬声器和已校准的麦克风。
记录了30+线性扫描150KHz-20KHz,长20ms,并以500记录了KS / s。
使用下面的Matlab / Octave脚本对齐和平均后的信号,在该脚本下可以看到结果信号。

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


X和XH并进行了上述计算,结果似乎是合理的。下面是H(紫色)和X/H(绿色)的归一化图。




该图已被截断为相关频率。

请让我知道我是否以错误的方式进行操作。域,我假设这将是一个简单的X/Hifft(X./H),但是到目前为止,我的所有尝试都未能获得任何合理的答案。可以在此处以mat7-binary格式找到频率向量wavwriteHfH

也许我只是累了,这里有一个简单的解决方案,但此刻我看不到它。任何帮助/建议都非常感谢。

评论

这个线程-dsp.stackexchange.com/questions/953/…-和这个dsp.stackexchange.com/questions/2705/…-对您可能有用。

是的,发现我的错误,谢谢吉姆。我只考虑信号的幅度,这导致了零相位时间信号。看来我现在得到了正确的结果,我将其添加为答案。

#1 楼

在查看吉姆·克莱(Jim Clay)在评论中提到的参考文献后找到了答案,谢谢吉姆。

我犯了一个错误,仅考虑幅度会导致零相位信号并且无法合理使用对于发射,至少在此设置中没有。

我最终使用的代码可以在下面看到。

该脚本遵循保留时域信号的命名约定


% 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
x conv h的频谱图如下所示:



尽管去卷积信号中有一些噪声,但这些对我来说似乎很合理。

下一个测试将看是否发射x deconv h除非扬声器无法发出那些频率,否则类似x_deconv_y

用测试结果更新

我们使用对数下扫来重述上述测量。这些结果似乎表明该方法可行。

验证测试包括发射x并期望取回X / H,即在所有频率下能量相等。
由于最差的输出频率约为比最好的要弱20dB,最高的输出电平会低得多。

发出的信号:



所记录信号的时间序列和频谱图如下所示:



评论


$ \ begingroup $
对此有任何更新吗?您发出的信号看起来如何?
$ \ endgroup $
–兰斯
2012年8月1日14:44

$ \ begingroup $
@Busk:感谢您的关注。由于该设备正在其他地方使用,因此我还没有机会对其进行测试。完成测试后,我将发布结果。
$ \ endgroup $
–雷神
2012年8月1日在16:04



$ \ begingroup $
@Busk:我们现在已经对其进行了测试,它似乎可以工作:-)。
$ \ endgroup $
–雷神
2012年10月31日8:33