信息处理 - 超声波扬声器校准和发射校准信号 - 吾爱随笔录

超声波扬声器校准和发射校准信号

信息处理 matlab 声音的频率响应

2022-01-06 13:36:36

我正在尝试校准超声波扬声器，目的是发出可预测的信号，但可惜我一直遇到麻烦，可能是因为我缺乏 DSP-fu。

一点背景

我希望能够播放尽可能接近校准录音的声音。据我了解这个理论，我需要找到扬声器的传递函数并对我想用它发出的信号进行反卷积。像这样的东西（在频域中）：

X -> H -> XH

发出X的信号H是扬声器的传递函数，XH是X次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/H我需要将其转换回时域之后，我认为这将是一个简单的ifft(X./H)and wavwrite，但到目前为止我的所有尝试都未能得到任何合理的答案。频率向量Hf,H和X可以在这里找到mat7-binary 格式。

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

1个回答

在查看了 Jim Clay 在评论中提到的参考资料后找到了答案，谢谢 Jim。

我犯了一个错误，即只考虑导致零相位信号的幅度并且不能明智地用于发射，至少在这个设置中没有。

我最终使用的代码如下所示。

该脚本遵循保持时域信号小写和频域信号大写的命名约定。

% 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 conv h 的频谱图 x deconv h 的频谱图