自适应滤波器做啸叫检测

%啸叫频点的检测
%啸叫位置不确定时的抑制方法
%初始条件设置
fs=44100;           %采样频率：44.1kHz
duration=1;         %音频时长1s
%保存的语音信号路径及名称
filename_howling_random='D:\\project_record\\speech_random_howling.wav';
filename_adapfilter_random='D:\\project_record\\speech_adapfilter_howling.wav';
speech_howling=audioread(filename_howling_random);
video_len=length(speech_howling);

%自适应滤波器需要参考信号，可以是标准的音频信号也可以是误差信号的线性变换
%本次输入的参考信号是误差信号
%因此最终的输出信号也是误差信号
%输出的误差才是滤波后的信号

[speech_howling,~]=audioread(filename_howling_random);
%信号截取的参数设定
%设定更新速度：
%每20ms更新（查看资料介绍，一般10ms~20ms更新速度较为适宜）
updRate=floor(20*fs/1000);         
 %设置帧长，单位ms
%40ms一帧，这样一帧内的音频信号可以近似看作是平稳的信号
t_zhen=40; 
fRate=fix(t_zhen*fs/1000);            
%求取样本点个数
n_samples=length(speech_howling);
%计算帧数
nFrames=fix(n_samples/updRate)-1; 
%初始化参数
n=1;
flag_zhen=0;
max_temp=0;
for t=1:nFrames
   if(n+fRate<=n_samples)
       zhen=speech_howling(n:n+fRate-1);
       Y_zhen=abs(fft(zhen));
       P_zhen=Y_zhen.^2;
       P_zhenaver=mean(P_zhen);
       [P_max,P_index]=max(P_zhen);
       if(10*log10(P_max)-10*log10(P_zhenaver)>10)
           if(max_temp==P_index)
             flag_zhen=flag_zhen+1;
           else
            flag_zhen=0;
            max_temp=P_index;
           end
       end    
       n_zhen=0:(fRate-1);
       N_zhen=fRate;
       f_zhen=n_zhen*fs/N_zhen;
       %每一帧画图，若想观察每一帧的频谱则取消注释
       %figure(t+1)
       %plot(f_zhen(1:fix(N_zhen/2)),Y_zhen(1:fix(N_zhen/2)));
       n=n+updRate;
   end
end
if(flag_zhen>=15)
    f_howling=(max_temp-1)*1000/t_zhen;
    disp('howling frequency is(Hz):')
    disp(f_howling);
end


%自适应滤波器的初始化 
%定义FIR滤波器阶数，一个频点需要二阶
M = 2;                                   
step = 0.1;                                %算法调节步长控制因子
y_out = zeros(video_len,1);              %滤波器输出
error_out = zeros(video_len,1);          %误差输出            
w_out = zeros(video_len,M);              %系数输出
t=(0:1:video_len-1)/video_len;%时间轴单位s

for i=1:video_len
    %数据输入
    if i == 1           %如果是第一次进入
        w = zeros(M,1); %初始化滤波器抽头系数
        x = zeros(M,1); %初始化信号向量
    end
    d = speech_howling(i);                     %输入新的期望信号
    x = [sin(2*pi*f_howling*(i-1)/fs)
         cos(2*pi*f_howling*(i-1)/fs)];                 %输入新的信号矢量
    %算法正体
    y = x' * w;                              %计算滤波器输出
    error = d - y;                           %计算误差
    w_forward = w + step * error * x;         %计算滤波器系数向量
    %变量更替
    w = w_forward;
    %滤波结果存储
    %其中y_out是啸叫信号
    %error_out是滤波后的音频信号
    %w_out是二阶滤波器参数
    y_out(i) = y;
    error_out(i) = error;
    w_out(i,:) = w';
end
 audiowrite(filename_adapfilter_random,error_out,fs);
sound(error_out,fs);
figure;
%=======时域图像======
t=(0:1:video_len-1)/video_len*duration;%时间轴单位s
subplot(2,1,1);
plot(t,error_out);
title('滤波后音频信号时域图像');
xlabel('时间：s');
grid on;
%=======频域图像======
n=0:video_len-1;        
N=video_len;           
Y1=fft(error_out,N);    
mag=abs(Y1);            
f=n*fs/N;               
%取1/2作图
subplot(2,1,2);
plot(f(1:fix(N/2)),mag(1:fix(N/2)));
title('滤波后音频信号频谱图');
xlabel('频率/Hz');      
ylabel('幅度');         
grid on;    
%画出FIR滤波器的零极点图