基于希尔伯特变换的移频方法用于抑制啸叫

最编程 2024-03-27 19:52:37

...

一、什么是啸叫（Howling）

在本地扩声系统中，麦克风拾取扬声器播放的音频后又传回给扬声器播放，当扩音的增益足够大，在某些频率就会产生自激振荡，从而产生啸叫。

图1 声音反馈回路，x(n)为麦克风采集的音频，y(n)为扬声器播放的音频

\frac{Y(\omega)}{S(\omega)} =\frac{G(\omega)}{1-G(\omega)F(\omega)}

根据奈奎斯特稳定判据，在某些频点满足以下条件将产生自激振荡：

G(\omega)F(\omega)>1

\measuredangle G(\omega)F(\omega)=2kπ

自激振荡导致该频点的信号幅度不断增大，最终形成啸叫。
啸叫抑制主要有相移/频移法、陷波法和自适应法，下面介绍基于希尔伯特变换的相移/频移法。

二、相移器

相移器的原理图如图2所示，其中粗线代表复数信号。输入信号是 $cos(ωt+ϕ)$ ，希尔伯特变换（Hilbert transformer）会将输入信号转换为 $e^{j(ωt+ϕ) } = cos(ωt+ϕ) + jsin(ωt+ϕ)$ 。将该信号与 $e^{jθ}$ （θ是期望的相位偏移）相乘，得到复数信号 $e^{j(ωt+ϕ+θ)}$ 。取该复数信号的实部就得到我们想要的输出： $y= cos(ωt+ϕ+θ)$ 。

图2 相移器原理

将希尔伯特复数输出表示为

I+jQ

，可以得到

y=Re\{(I+jQ)e^{jθ} \}=Re\{(I+jQ)(cosθ+jsinθ)\}

y=Icosθ-Qsinθ \qquad \qquad (公式1)

从上式可得到相移器的实现（见图3）。注意图3中没有出现“j”，这意味着处理的信号都是实信号！

图3 相移器的实现

下面用Matlab代码展示图3所示的相移器的实现。31抽头的希尔伯特变换器的实现如下面代码所示，对系数的理论值（第4节将介绍如何计算理论值）乘以Hamming窗后得到系数b1，b2是15个采样点延迟（对应31抽头希尔伯特变换器中心抽头的延迟）。

    % 31-tap Hilbert transformer
    b= 2/pi * [-1/15 0 -1/13 0 -1/11 0 -1/9 0 -1/7 0 -1/5 0 -1/3 0 -1 0 1 0 … 
                1/3 0 1/5 0 1/7 0 1/9 0 1/11 0 1/13 0 1/15];
    b1= b.*hamming(31)';                          % window the coefficients
    b2= [0 0 0 0  0 0 0 0  0 0 0 0  0 0 0 1];     % delay of 15 (center tap of HT)

接下来创建一个6Hz的正弦输入信号 $x$ ，其采样率为100Hz：

    fs= 100;                    % Hz sample frequency
    Ts= 1/fs;
    N= 128;
    n= 0:N-1;
    f0= 6;                       % Hz freq of input signal
    x= cos(2*pi*f0*n*Ts);        % input signal

现在对 $x$ 执行希尔伯特变换：

    I= filter(b2,1,x);        % I= center tap of HT
    Q= filter(b1,1,x);        % Q= output of HT

最后根据公式1得到相移 $–π/3$ 的输出：

    theta= -pi/3;                          % rad phase shift with respect to I
    y= I*cos(theta) - Q*sin(theta);        % phase shifter output

在图3中我们假设希尔伯特变换不会引起相移，但事实上在实现中会引入15个采样点的延迟，换句话说我们是在希尔伯特转换输出的 $I$ 的基础上增加 $–π/3$ 的相移得到输出 $y$ 。图4打印了 $y$ 和 $I$ 的32个样点，可以看到 $y$ 和 $I$ 相位相差 $π/3$ 。

图4 相移器的输入和输出，f0=6Hz，fs=100Hz，ϕ= - π/3，蓝线=I，绿线=y

三、频移器

频移器的实现如图5所示，它和相移器类似，只是我们用 $dω*t$ 替换公式1中的 $θ$ ，其中 $dω$ 为期望的频率偏移。

图5 频移器的实现

下面展示将中心频率为12Hz的输入信号频移+3Hz的Matlab代码。首先将相移器例子中展示的希尔伯特变换器的Hamming窗替换为Blackman窗。与Hamming窗相比，Blackman窗有更精确的通带相应，代价则是牺牲了低频的分辨率。接下来将参数量化到小数点后12bits。

    % 31-tap Hilbert transformer
    b= 2/pi * [-1/15 0 -1/13 0 -1/11 0 -1/9 0 -1/7 0 -1/5 0 -1/3 0 -1 0 1 0 … 
               1/3 0 1/5 0 1/7 0 1/9 0 1/11 0 1/13 0 1/15];
    b1= b.*blackman(31)';                         % window the coefficients
    b1= round(b1*2^12)/2^12;                      % quantize coefficients
    b2= [0 0 0 0  0 0 0 0  0 0 0 0  0 0 0 1];     % delay of 15 (center tap of HT)

生成输入信号 $x$ ，并对它执行希尔伯特变换。这次的输入信号不用正弦信号，而是选择有矩形频谱的脉冲调制信号。

    fs= 100;                       % Hz sample frequency
    Ts= 1/fs;
    N= 2048;
    n= 0:N-1;
    fc= 12;                        % Hz carrier frequency
    bw= 2;                         % Hz -3 dB bandwidth of modulated pulse
    x= .5*modpulse(fc,bw,N,fs);    % modulated pulse with approx rect spectrum
    % Apply modulated pulse to Hilbert transformer
    I= filter(b2,1,x);             % I= center tap of HT filter
    Q= filter(b1,1,x);             % Q= output of HT filter

接着通过NCO（数字振荡器）生成正弦信号与余弦信号，并利用它们实现频移。

    df= 3;                              % Hz desired frequency shift
    u= mod(df*n*Ts,1);                  % phase accumulator
    theta= 2*pi*u;                      % phase
    y= I.*cos(theta) –Q.*sin(theta);    % final output

现在我们已经可以得到实数输出 $y$ 的频谱了，但如果我们将图5中所有的信号放在一起展示会看到更多有趣的内容。图5有实数信号 $x$ 、 $y$ 和三个复数信号：希尔伯特变换的输出、NCO的输出、这两者相乘的输出。这些复数信号可以通过以下代码得到：

    xc= I + j*Q;            % complex Hilbert transformer output
    nco= exp(j*theta);      % complex nco output
    yc= xc.*nco;            % complex product

获取这些信号的频谱：

    X= fft(x,N);          XdB= 20*log10(abs(x));
    XC= fft(xc,N);        XCdB= 20*log10(abs(XC));
    YC= fft(yc,N);        YCdB= 20*log10(abs(YC));
    NCO= psd(nco,N,fs);   NCOdB= 10*log10(abs(NCO));
    Y= fft(y,N);          YdB= 20*log10(abs(Y));

这些频谱在图6中展示。其中图6a是实数输入 $x$ 的频谱。图6b是希尔伯特变换的输出 $xc$ ，可看作解析信号（analytic signal，没有负频率分量的复数信号），它的负频率分量几乎可以忽略。图6c是NCO信号 $e^{jθ}$ ，同样是个解析信号，只在+3Hz处有值。图6d是 $xc$ 和 $e^{jθ}$ 相乘的结果 $yc$ ，可以看到两者相乘的结果是 $xc$ 的中心频率从12Hz移到15Hz。最后图6e是实数输出 $y$ ，中心频率为15Hz。
假如我们只是对两个实数信号相乘，而不是对两个复数信号相乘，那得到的会是中心频率为9Hz和15Hz的两个信号叠加的频谱，这并非我们期望的结果。

图6 a.实数输入X b.希尔伯特输出XC c.NCO输出 d.相乘的结果YC e.最终输出Y

四、希尔伯特变换

下面介绍希尔伯特变换的基本原理和实现。
一个理想的离散希尔伯特变换的频率响应为： $H(z)=-j,\qquad 0<f<f_{s}/2$ $\qquad \qquad \qquad \qquad=j,\qquad -f_{s}/2<f<0 \qquad (公式2)$
该响应可以通过一个滤波器近似得到，该滤波器的脉冲响应为：
$h[n]=\frac{2}{πn}, \qquad n为奇数$ $\qquad \qquad \qquad=0, \qquad n为偶数 \qquad (公式3)$
该脉冲响应是非因果的，但可以通过截断来近似，令 $n=-N/2:N/2-1$ ，其中 $N$ 为偶数，且 $N+1$ 为抽头数。抽头系数可以乘上一个窗函数。图7、图8显示7抽头希尔伯特变换器，其中图7抽头系数序号与公式3保持一致，而图8重排了序号，让抽头系数从 $b_0$ 开始。 $I$ 通道在滤波器延迟网络的正中间，对应着理想希尔伯特滤波器的 $t=0$ 时刻。
Q通道的响应 $-j=e^{-jπ/2}$ 对应 $π/2$ 的相位滞后。因此对于 $I=cos(ωt)$ ，相位滞后 $π/2$ ，得到 $Q$ 为 $sin(ωt)$ 。复数信号 $I + jQ = cos(ωt) +jsin(ωt) = e^{jωt}$ 。与实数输入信号不同，该信号没有负频率分量，是一个解析信号。

图7 抽头系数序号与公式3保持一致

图8 抽头系数从b0开始

在前面的章节我们实现了一个31抽头的希尔伯特变换器，代码如下：

     % 31-tap Hilbert transformer
    b= 2/pi * [-1/15 0 -1/13 0 -1/11 0 -1/9 0 -1/7 0 -1/5 0 -1/3 0 -1 0 1 0 … 
               1/3 0 1/5 0 1/7 0 1/9 0 1/11 0 1/13 0 1/15];
    b1= b.*blackman(31)';                         % window the coefficients
    b1= round(b1*2^12)/2^12;                      % quantize coefficients
    b2= [0 0 0 0  0 0 0 0  0 0 0 0  0 0 0 1];     % delay of 15 (center tap of HT)

这些滤波器系数如图9所示。下面计算频率响应。为了得到一个可实现的滤波器，我们对输入信号延迟了15个采样点得到 $I$ 通道，因此我们要计算的频率响应也是延迟了15个采样点的输入信号的响应。下面Matlab代码计算得到的频率响应如图10所示，可以看到在低频和频率接近 $f_s/2$ 的地方，频率响应与公式2的理想希尔伯特变换的频率响应有较大偏差，因此我们在执行希尔伯特变换的时候，要确保输入信号的频率落在图10中频率响应的平坦部分。

    N= 1024;
    k= -N/2:N/2-1;
    f= k*fs/N;
    h_Q= fft(b1,N);            % h_Q = response of b1
    h_I= fft(b2,N);            % h_I = response at center tap (delay of 15 samples)
    h= h_Q./h_I;
    h= fftshift(h);            %shift dc to center (swap left and right halves of h)
    plot(f,imag(h))

图9 31抽头的希尔伯特变换器系数

图10 31抽头的希尔伯特变换器的频率响应，fs=100Hz

参考文章
Phase or Frequency Shifter Using a Hilbert Transformer
啸叫抑制之陷波法

上一篇：相位和相位差

下一篇： m 基于 FPGA 的 DQPSK 调制解调通信系统 verilog 实现，包括测试平台，不包括载波同步