音頻進階學習十九——逆系統（簡單進行回聲消除）

文章目錄

前言
一、可逆系統
- 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;