音频进阶学习十九——逆系统（简单进行回声消除）

文章目录

前言
一、可逆系统
- 1.定义
- 2.解卷积
- 3.逆系统恢复原始信号过程
- 4.逆系统与原系统的零极点关系
二、使用逆系统去除回声
- 获取原信号的频谱
- 原系统和逆系统幅频响应和相频响应
- 使用逆系统恢复原始信号
- 整体代码如下
总结

前言

在上一篇音频进阶学习十八——幅频响应相同系统、全通系统、最小相位系统文章中，我们在最小相位系统中提到了逆系统：一个稳定因果的LTI系统，同时会有一个稳定因果的逆系统。

本文中将会介绍逆系统的含义，通过对于零极点的分析，来对比逆系统和原系统的关联。最后，将会使用逆系统来对于给定的近端信号和参考信号来进行回声消除。

|版本声明：山河君，未经博主允许，禁止转载

一、可逆系统

1.定义

如果一个系统的输入是 $x [n]$ ，经过系统 $H (z)$ 后系统的输出是 $y [n]$ ，那么将 $y [n]$ 通过系统 $G (z)$ 后输出的是 $x [n]$ ，此时系统 $G (z)$ 是系统 $H (z)$ 的逆系统。

而 $G (z)$ 和 $H (z)$ 的关系：

频域关系
$G (z) H (z) = 1 => g [n] = h [n]$
时域关系
对于系统 $H$ 的脉冲响应有： $y [n] = x [n] * h [n]$
对于系统 $G$ 的脉冲响应有： $x^{'} [n] = y [n] * g [n] = (x [n] * h [n]) * g [n] = x [n] * (h [n] * g [n])$
如果想要 $x^{'} [n] = x [n]$ ，那么 $h[n]*g[n]=\delta n$

2.解卷积

对于给定的输入信号 $x [n]$ 和输出信号 $y [n]$ ，我们知道通过系统有：
$X(f)=F\{x[n]\},Y(f)=F\{x[n]\},H(f)=\frac{X(f)}{Y(f)}$
通过逆变换， $h[n]=F^{-1}\{H(f)\}$ 得到冲激响应的这一过程，叫做解卷积

3.逆系统恢复原始信号过程

首先结合具体场景，假设发射了一个信号 $x [n]$ ，通过多条途径例如墙体、窗户的反射，导致接收的信号是 $y [n]$ ，我们想通过 $y [n]$ 恢复为原始的 $x [n]$ ，如下图：
在这里插入图片描述

发射约定好的信号 $T [n]$
接收 $T^{'} [n] = T [n] * h [n]$
通过解卷积得到 $h [n]$
构建逆系统 $g [n]$
通过逆系统获得原始信号 $x [n] = y [n] * G [n]$

值得注意的是，在实际场景中，过程会比较复杂，例如不存在发射约定好的信号，或者场景变换导致冲激响应的变化等等，这就需要一个自适应滤波器不断的估计回声信号，当然这是以后的内容。

4.逆系统与原系统的零极点关系

对于Z变换的频率响应，它的逆系统
$H(z)=\frac{\prod_{m=1}^M(1-c_mz^{-1})}{\prod_{n=1}^N(1-d_nz^{-1})},G(z)=\frac{\prod_{n=1}^N(1-d_nz^{-1})}{\prod_{m=1}^M(1-c_mz^{-1})}$
也就是原系统 $H (z)$ 和逆系统 $G (z)$ 的零极点是互换的，那么此时如果逆系统 $G (z)$ 是因果稳定的，那么逆系统 $G (z)$ 的零点和极点都在单位圆内，而此时原系统的 $H (z)$ 零点和极点同样都在单位圆内，此时都为最小相位系统。

这也是上一篇文章中最小相位系统提到的特性。

二、使用逆系统去除回声

现在使用matlab来设计一个逆系统，对于已知的输入信号和输出信号进行转换。

获取原信号的频谱

需要注意的是，这里的输入信号和输出信号是笔者自己通过pcm叠加转换的，所以频域上可能看不出特别大的区别，但是实际可以明显听到回声。

核心代码

%输入信号和输出信号频域图
H_in=fft(input_signal);
H_out=fft(output_signal);
R_in = abs(H_in);
R_out = abs(H_out);
freq=linspace(0,fs/2, floor(min_len/2));figure;
subplot(2,1,1);
plot(freq, R_in(1:floor(min_len/2)));
title("输入信号频域图");xlabel("频率H(z)");ylabel("幅度");grid on;
subplot(2,1,2);
plot(freq, R_out(1:floor(min_len/2)));
title("输出信号频域图");xlabel("频率H(z)");ylabel("幅度");grid on;

得到结果：
在这里插入图片描述

原系统和逆系统幅频响应和相频响应

通过图像幅频响应图像和相频响应的图像，可以很清楚的看到原系统和逆系统的幅频响应和相频响应是倒过来的。

同时值得注意的是：

某些频率上幅值接近零（即原系统在这些频率上有很强的衰减），那么其倒数 $H_{inv}$ 可能会变得非常大，导致放大效应。
在计算 $H_{inv}$ 浮点数精度问题可能会导致不稳定。
幅度部分可能被正常恢复，但相位的误差会导致合成信号在时域上发生相干叠加，造成某些时刻的峰值更大，最终影响整体信号能量。

而这些问题需要结合窗口平滑来进行解决，这里是通过增加了一些阈值来减少极端值的出现。

代码如下：

H_c = H_in ./ (H_out + eps);
%H_inv= 1 ./ (H_c + eps);
%H_inv(abs(H_c) < 1e-3) = 0;
H_inv = 1 ./ (H_c + 0.01);
H_inv(abs(H_c) < 1e-2) = 0;
R_c=abs(H_c);
R_inv=abs(H_inv);figure;
subplot(2,2,1);
plot(freq, 20*log10(R_c(1:floor(min_len/2))));
title('原系统幅频响应');xlabel('频率 (Hz)'); ylabel('幅度 (dB)'); grid on;
subplot(2,2,2);
plot(freq, rad2deg(angle(H_c(1:floor(min_len/2)))));
title('原系统相频响应');xlabel('频率 (Hz)'); ylabel('相位(度)'); grid on;
subplot(2,2,3);
plot(freq, 20*log10(R_inv(1:floor(min_len/2))));
title('逆系统幅频响应');xlabel('频率 (Hz)'); ylabel('幅度 (dB)'); grid on;
subplot(2,2,4);
plot(freq, rad2deg(angle(H_inv(1:floor(min_len/2)))));
title('逆系统相频响应');xlabel('频率 (Hz)'); ylabel('相位(度)'); grid on;

得到
在这里插入图片描述

使用逆系统恢复原始信号

最后恢复的信号还是变大了，但是此时已经听不见回声了，这里在恢复时已经使用整体能量归一化来尽量减少能量放大的影响。

代码如下：

%还原后的信号频谱
restored_signal = real(ifft(fft(output_signal) .* H_inv));
restored_signal = restored_signal / norm(restored_signal) * norm(input_signal);restored_signal = int16(restored_signal);
fwrite(turn_file, restored_signal, 'int16');
fclose(turn_file);
H_restored = fft(restored_signal);
R_res = abs(H_restored);
figure;
subplot(1,1,1);
plot(freq, R_res(1:floor(min_len/2)));
title('还原信号的频域图');
xlabel('频率 (Hz)'); ylabel('幅度'); grid on;

得到
在这里插入图片描述

整体代码如下

fs=16000;
input_file=fopen("input.pcm",'rb');
output_file=fopen("output.pcm",'rb');
turn_file=fopen("turn.pcm",'wb+');input_signal=fread(input_file, inf, 'int16');
output_signal=fread(output_file, inf, 'int16');
fclose(input_file);
fclose(output_file);min_len=min(length(input_signal), length(output_signal));
input_signal=input_signal(1:min_len);
output_signal=output_signal(1:min_len);%输入信号和输出信号频域图
H_in=fft(input_signal);
H_out=fft(output_signal);
R_in = abs(H_in);
R_out = abs(H_out);
freq=linspace(0,fs/2, floor(min_len/2));figure;
subplot(2,1,1);
plot(freq, R_in(1:floor(min_len/2)));
title("输入信号频域图");xlabel("频率H(z)");ylabel("幅度");grid on;
subplot(2,1,2);
plot(freq, R_out(1:floor(min_len/2)));
title("输入信号相频响应");xlabel("频率H(z)");ylabel("幅度");grid on;%原系统和逆系统幅频响应和相频响应
H_c = H_in ./ (H_out + eps);
%H_inv= 1 ./ (H_c + eps);
%H_inv(abs(H_c) < 1e-3) = 0;
H_inv = 1 ./ (H_c + 0.01);
H_inv(abs(H_c) < 1e-2) = 0;
R_c=abs(H_c);
R_inv=abs(H_inv);figure;
subplot(2,2,1);
plot(freq, 20*log10(R_c(1:floor(min_len/2))));
title('原系统幅频响应');xlabel('频率 (Hz)'); ylabel('幅度 (dB)'); grid on;
subplot(2,2,2);
plot(freq, rad2deg(angle(H_c(1:floor(min_len/2)))));
title('原系统相频响应');xlabel('频率 (Hz)'); ylabel('相位(度)'); grid on;
subplot(2,2,3);
plot(freq, 20*log10(R_inv(1:floor(min_len/2))));
title('逆系统幅频响应');xlabel('频率 (Hz)'); ylabel('幅度 (dB)'); grid on;
subplot(2,2,4);
plot(freq, rad2deg(angle(H_inv(1:floor(min_len/2)))));
title('逆系统相频响应');xlabel('频率 (Hz)'); ylabel('相位(度)'); grid on;%还原后的信号频谱
restored_signal = real(ifft(fft(output_signal) .* H_inv));
restored_signal = restored_signal / norm(restored_signal) * norm(input_signal);restored_signal = int16(restored_signal);
fwrite(turn_file, restored_signal, 'int16');
fclose(turn_file);
H_restored = fft(restored_signal);
R_res = abs(H_restored);
figure;
subplot(1,1,1);
plot(freq, R_res(1:floor(min_len/2)));
title('还原信号的频域图');
xlabel('频率 (Hz)'); ylabel('幅度'); grid on;