一、什么是啸叫（Howling）

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

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

二、相移器

相移器的原理图如图2所示，其中粗线代表复数信号。输入信号是，希尔伯特变换（Hilbert transformer）会将输入信号转换为。将该信号与（θ是期望的相位偏移）相乘，得到复数信号。取该复数信号的实部就得到我们想要的输出：。

图2 相移器原理

图3 相移器的实现

    % 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的正弦输入信号，其采样率为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

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

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

最后根据公式1得到相移的输出：

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

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

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

三、频移器

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

图5 频移器的实现

    % 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)

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

    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

现在我们已经可以得到实数输出的频谱了，但如果我们将图5中所有的信号放在一起展示会看到更多有趣的内容。图5有实数信号、和三个复数信号：希尔伯特变换的输出、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是实数输入的频谱。图6b是希尔伯特变换的输出，可看作解析信号（analytic signal，没有负频率分量的复数信号），它的负频率分量几乎可以忽略。图6c是NCO信号，同样是个解析信号，只在+3Hz处有值。图6d是和相乘的结果，可以看到两者相乘的结果是的中心频率从12Hz移到15Hz。最后图6e是实数输出，中心频率为15Hz。
假如我们只是对两个实数信号相乘，而不是对两个复数信号相乘，那得到的会是中心频率为9Hz和15Hz的两个信号叠加的频谱，这并非我们期望的结果。

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

四、希尔伯特变换

下面介绍希尔伯特变换的基本原理和实现。
一个理想的离散希尔伯特变换的频率响应为：
该响应可以通过一个滤波器近似得到，该滤波器的脉冲响应为：

该脉冲响应是非因果的，但可以通过截断来近似，令，其中为偶数，且为抽头数。抽头系数可以乘上一个窗函数。图7、图8显示7抽头希尔伯特变换器，其中图7抽头系数序号与公式3保持一致，而图8重排了序号，让抽头系数从开始。通道在滤波器延迟网络的正中间，对应着理想希尔伯特滤波器的时刻。
Q通道的响应对应的相位滞后。因此对于，相位滞后，得到为。复数信号。与实数输入信号不同，该信号没有负频率分量，是一个解析信号。

图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个采样点得到通道，因此我们要计算的频率响应也是延迟了15个采样点的输入信号的响应。下面Matlab代码计算得到的频率响应如图10所示，可以看到在低频和频率接近的地方，频率响应与公式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
啸叫抑制之陷波法