啸叫抑制（AFS）从算法仿真到工程源码实现-第七节-自适应滤波法

一、概述

        自适应滤波器做啸叫抑制，类似回声消除的方法，刚开始受一些调研资料的影响，由于信号的相关性，自适应滤波器会完全消除语音。后来随着研究的深入，语音帧之间的相关性没有想象中的强烈，是当前主流的方法。直接使用NLMS算法在时域就可以取得较好的效果。后面进行了工程化，使用频域分块自适应滤波器，取得了理想的效果。关于自适应滤波器进行啸叫抑制，参考数据必须是扬声器播放之前的数据，也就是自动增益控制算法之后的数据。

二、算法仿真

2.1 算法流程图

2.2 时域python仿真

#!/usr/bin/python
from __future__ import division
# import

import numpy as np
import matplotlib.pyplot as plt
import soundfile as sf


def simulate_howling_env():
    input_file = "test/LDC93S6A.wav"
    howling_file = "test/remove_howling.wav"

    # load clean speech file
    x, Srate = sf.read(input_file)

    # pre design a room impulse response
    rir = np.loadtxt('test/path.txt', delimiter='\t')

    # G : gain from mic to speaker
    G = 0.2

    # ====== set parameters ========
    interval = 0.02  # frame interval = 0.02s
    Slen = int(np.floor(interval * Srate))
    if Slen % 2 == 1:
        Slen = Slen + 1
    PERC = 50  # window overlap in percent of frame size
    len1 = int(np.floor(Slen * PERC / 100))
    len2 = int(Slen - len1)
    nFFT = 2 * Slen

    N = min(2000, len(rir))  # limit room impulse response length
    x2 = np.zeros(N)  # buffer N samples of speaker output to generate acoustic feedback
    y = np.zeros(len(x))  # save speaker output to y
    y1 = 0.0  # init as 0

    w = np.zeros(N)
    u = np.zeros(N)

    mu = 0.01
    for i in range(len(x)):
        x1 = x[i] + y1

        e = x1 - np.dot(u, w)
        w = w + mu * e * u / (np.dot(u, u) + 1e-6)
        x1 = e

        y[i] = G * x1
        y[i] = min(2, y[i])  # amplitude clipping
        y[i] = max(-2, y[i])
        x2[1:] = x2[:N - 1]
        x2[0] = y[i]

        u[1:] = u[:N - 1]
        u[0] = y[i]

        y1 = np.dot(x2, rir[:N])

    sf.write(howling_file, y, Srate)
    plt.figure()
    plt.subplot(2, 1, 1)
    plt.plot(y)
    plt.xlim((0, len(y)))
    plt.subplot(2, 1, 2)
    plt.specgram(y, NFFT=nFFT, Fs=Srate, noverlap=len2, cmap='jet')
    plt.ylim((0, 5000))
    plt.ylabel("Frquency (Hz)")
    plt.xlabel("Time (s)")
    plt.show()


simulate_howling_env()

2.3 时域c仿真

#include <sys/types.h>
#include <sys/stat.h>
#include <stdlib.h>
#include <stdio.h>
#include <fcntl.h>
#include <unistd.h>
#include <string.h>

#define OUTPUT_PATH_NAME "/tmp/node1_to_node2.tmp"
#define INPUT_PATH_NAME "/tmp/node2_to_node1.tmp"

int main() {

  if (access(INPUT_PATH_NAME, F_OK) != -1) {
    printf("File exists.\n");
  } else {
    mkfifo(INPUT_PATH_NAME, 0666);
  }
  int fin = open(INPUT_PATH_NAME, O_RDONLY);
  int fout = open(OUTPUT_PATH_NAME, O_WRONLY);
  int frames = 128;
  int size = frames * 2;
  char *buffer = (char *) malloc(size);
  
  int chunk_num = 5;
  float mu = 0.01;

  float *u = (float*) malloc(sizeof(float) * frames * chunk_num);
  float *w = (float*) malloc(sizeof(float) * frames * chunk_num);
  
  for (int i = 0; i < frames * chunk_num; i++) {
    u[i] = 0;
	w[i] = 0;
  }

  float *current_out_data = (float*) malloc(sizeof(float) * frames);
  float *current_data = (float*) malloc(sizeof(float) * frames);
  float *last_data = (float*) malloc(sizeof(float) * frames);
  
  short *out_data = (short*) malloc(sizeof(short) * frames);
  
  for (int i = 0; i < frames; i++) {
    current_out_data[i] = 0;
    current_data[i] = 0;
    last_data[i] = 0;
	out_data[i] = 0;
  }

  while (1) {
    int rc = read(fin, buffer, size);
    if (rc <= 0) {
      fprintf(stderr, "end of file on input/n");
      break;
    } else if (rc != size) {
      fprintf(stderr,
        "short read: read %d bytes/n", rc);
    } else {
      printf("read len: %d\n", rc);
    }
	
	short *input_data = (short*)buffer;
	//printf("current_data:");
	for (int i = 0; i < frames; i++) {
	  current_data[i] = (float)(input_data[i]) / 32768.0;
	  //printf(" %f", current_data[i]);
	}
	//printf("\n");
	
	//printf("current_out_data:");
	for (int i = 0; i < frames; i++) {
	  memmove(u+1, u, sizeof(float)*(frames*chunk_num-1));
	  u[0] = last_data[i];
	  
	  float uw = 0;
	  float uu = 1e-6;
	  for (int k = 0; k < frames*chunk_num; k++) {
	    uw += u[k]*w[k];
		uu += u[k]*u[k];
	  }
	  
	  current_out_data[i] = current_data[i] - uw;
	  //printf(" (%f %f %f)", current_out_data[i], uw, uu);
	  
	  for (int k = 0; k < frames*chunk_num; k++) {
	    w[k] = w[k] + mu * current_out_data[i] * u[k] / uu;
	  }
	}
	//printf("\n");
	
	//printf("out_data:");
	for (int i = 0; i < frames; i++) {
	  if (current_out_data[i] > 1.0) {
	    current_out_data[i] = 1.0;
	  } else if (current_out_data[i] < -1.0) {
	    current_out_data[i] = -1.0;
	  }
	  out_data[i] = (short)((current_out_data[i]) * 32767.0);
	  //printf(" %d", out_data[i]);
	}
	//printf("\n");
	
	for (int i = 0; i < frames; i++) {
	  last_data[i] = current_data[i];
	}

    rc = write(fout, (char*)out_data, size);
    if (rc <= 0) {
      break;
    } else {
      printf("write len: %d\n", rc);
    }
  }
  close(fin);
  close(fout);
  free(buffer);
  
  free(u);
  free(w);
  
  free(current_out_data);
  free(current_data);
  free(last_data);
  
  free(out_data);
  
  unlink(INPUT_PATH_NAME);
  return 0;
}

2.4 仿真结果

三、总结

        本节我们使用自适应滤波法在时域进行去啸叫仿真，下节我们进行频域的实现，完成整个系统的搭建。