找到滤波器的“通道”或响应的最小均方解由以下 MATLAB/Octave 代码提供,使用滤波器的输入作为 tx,滤波器的输出作为 rx。有关其工作原理的更多详细信息,请参阅这篇文章:补偿音频信号中的扬声器频率响应:
function coeff = channel(tx,rx,ntaps)
% Determines channel coefficients using the Wiener-Hopf equations (LMS Solution)
% TX = Transmitted (channel input) waveform, row vector, length must be >> ntaps
% RX = Received (ch output) waveform, row vector, length must be >> ntaps
% NTAPS = Number of taps for channel coefficients
% Dan Boschen 1/13/2020
tx= tx(:)'; % force row vector
rx= rx(:)'; % force row vector
depth = min(length(rx),length(tx));
A=convmtx(rx(1:depth).',ntaps);
R=A'*A; % autocorrelation matrix
X=[tx(1:depth) zeros(1,ntaps-1)].';
ro=A'*X; % cross correlation vector
coeff=(inv(R)*ro); %solution
end
OP 使用梳状滤波器的情况是最具挑战性的情况之一,因为它取决于解决方案的每个频率处的信号能量(这就是为什么线性均衡器,如果你交换 rx 和 tx,这个函数正在执行,不执行在频率选择信道中表现良好,并导致在零位置处的噪声增强)。在使用 MATLAB 或 Octave 确定的测试滤波器的频率响应下方,显示与此类梳状滤波器相关的多个频率零点:
b=[1,zeros(1,39),-1];
freqz(b,1,2^14) % 2^14 samples to show nulls

MATLAB 或 Octave 脚本来演示上述函数的使用并确定输出和输入之间的延迟:
%% Filter with OP's example
b=[1,zeros(1,39),-1]; % numerator coefficients
a = 1; % denominator coefficients
%% Generate signal using OP's code
Fs=1000;
t=0:1/Fs:(5-1/Fs);
wi=blackman(length(t))';
rn=+rand(1,length(t))*.2;
x1=sin(2*pi*13*t).*wi +rn;
x2=sin(2*pi*25*t).*wi +rn;
x3=sin(2*pi*2*t).*wi +rn;
x=[x1,x2,x3];
% Filter
y=filter(b,a,x);
%% Test filter estimation
cf=channel(x,y,61);
%compare original and estimated channel
subplot(2,1,1)
stem(b)
title("Actual Channel Response")
xlabel("Sample Number")
subplot(2,1,2)
stem(cf)
title("Estimated Channel Response")
xlabel("Sample Number")

在调用函数通道时,我们可以准确地使用 41 次抽头,这样就可以正确解决解决方案;然而,最佳实践是从更长的滤波器长度开始,评估结果,然后相应地减少抽头。在噪声条件下的实际实践中,使用比必要更多的抽头会导致噪声增强,因此最好使用捕获主要抽头权重所需的最小抽头来完成最终解决方案。
使用 MATLAB 和 Octaves grpdelay 命令观察 groupdelay 图的问题,即无法解决没有信号通过滤波器的延迟问题(任何事情都很难确定在这些频率之一的单个音调的延迟,即被滤波器归零!)但它能够准确地确定信号能量存在的延迟。同样,波形本身必须在我们寻求解决方案的所有频率上都有能量。OP 的测试波形的频谱密度在所有频率上充分分布以适合此目的(这就是伪随机波形产生良好“通道探测”模式的原因)。
该图用于与 OP 的图进行比较,显示该滤波器的群延迟为 20 个样本。

这是另一个测试用例,展示了它如何与更合理的通道(无深空值)一起使用,其系数和频率响应低于此:
b = [0.2 .4 -.3 .4 .3 .1];

该解决方案在实际和估计之间无法区分,因此我添加了噪声以使其更有趣,使用上面代码中的相同 x 和 y:
% add noise
noise = 0.351*randn(1,length(y));
yn = y + noise;
snr = 20*log10(std(yn)./std(noise));
%% Test filter estimation
cf=channel(x,yn,10);
