实验二利用MATLAB工具箱对混杂噪声的音频信号进行滤波

实验二利用MATLAB工具箱对混杂噪声的音频信号进行滤波

一、小提琴音频的加噪去噪处理

要求：选择子作业１中的音频信号，自行给定滤波器的系统函数，分别采用时域线性卷积和差分方程两种方法对音频信号进行滤波处理，比较滤波前后信号的波形和回放的效果。
最终整体效果如下图：

1、音频信号的构建

①音乐信号的产生：

  采用Adobe Audition提取出一首小提琴音乐的一个音符的音频信号，如下图所示：

上下信号波形分别为左右声道。如何判断为一个音符的信号呢？通过Adobe Audition的频谱分析进行判断，语谱图如下图所示：

语谱图是将信号的频谱和时间结合，横轴为时间，纵轴为频率，颜色的深浅表征信号的幅度。通过在一段音乐信号的语谱图中，找到在一定时间段内各频率分量稳定的区间，则可判断其为一个音符。通过上图的各频率成分的颜色可得知，在该时间段下，在某几段频率区间中，颜色波形十分稳定。
  通过MATLAB对该音乐信号进行频谱分析可得：（左图为时域信号，右图为频域信号）

将上图与前一幅图对比可知，频域图显示能量高的地方，其对应语谱图上的颜色更鲜艳，则说明其表征了音频信号的能量分布。

②噪声信号的产生：

  因为本实验是要对信号作滤波处理，所以我在该音乐信号中通过Adobe Audition软件添加了主要能量分布在低频的高斯白噪声，其时域波形和语谱图如下图所示：

由上图可见，该噪声能量的大部分分布在400Hz以下，主要能量分布在100Hz以下，因此后续的滤波器设计选择使用低通滤波器。
  通过MATLAB对该噪声信号进行频谱分析可得：（左图为时域信号，右图为频域信号）

可见该噪声的主要能量分布在200Hz以下。

③混杂信号的产生：

  通过Adobe Audition软件的多轨功能将音乐信号和噪声信号进行叠加，因此该噪声为加性噪声，叠加后的时域图和语谱图如下图所示：

由上图可知，在200Hz的频率成分以下叠加了主要的噪声。
  通过MATLAB对该混杂信号进行频谱分析可得：（左图为时域信号，右图为频域信号）

2、滤波器的设计

①MATLAB滤波器设计工具：

  在MATLAB的应用程序中选择Filter Design & Analysis如下图所示：

则得到如下界面：

以下为对上图序号框区域的介绍：

1. 选择高通/低通/带通/带阻等滤波器类型以及滤波器实现方法（如IIR或FIR）；
2. 选择滤波器阶数，阶数越高，陡峭性越好；
3. 分别设置幅频特性坐标的单位，信号采样率，通带截止频率和阻带截止频率；
4. 该区域可以通过在Analysis菜单栏进行选择来显示该滤波器的不同参数指标或者数据分析；
5. 该区域为所设计滤波器的简要信息：结构，阶数，稳定性和滤波器来源。

②滤波器参数

  上图为我所设计的滤波器，为高通滤波器，IIR结构的巴特沃斯滤波器，阶数为由软件自动设置为符合条件的最小阶数45阶，原信号的采样率8000，通带截止频率为250Hz，阻带截止频率为200Hz。

在上图所示的选项卡中可以看到滤波器系数，这个系数为系统函数$H(Z)$分子和分母在$Z^{-n}$前面的系数，它是由23个分式相乘而成。
  该工具生成的代码如下所示：

    function Hd = my_filter_code2
    %MY_FILTER_CODE2 Returns a discrete-time filter object.
    
    % MATLAB Code
    % Generated by MATLAB(R) 9.0 and the Signal Processing Toolbox 7.2.
    % Generated on: 19-Oct-2022 14:56:51
    
    % Butterworth Highpass filter designed using FDESIGN.HIGHPASS.
    
    % All frequency values are in Hz.
    Fs = 8000;  % Sampling Frequency
    
    Fstop = 200; % Stopband Frequency
    Fpass = 250; % Passband Frequency
    Astop = 80;  % Stopband Attenuation (dB)
    Apass = 1;   % Passband Ripple (dB)
    match = 'stopband';  % Band to match exactly
    
    % Construct an FDESIGN object and call its BUTTER method.
    h  = fdesign.highpass(Fstop, Fpass, Astop, Apass, Fs);
    Hd = design(h, 'butter', 'MatchExactly', match);
    
    % [EOF]

3、利用该滤波器对混杂音频信号进行滤波

  主要代码如下：

    % ==========MATLAB工具箱生成的滤波器========== %
    H = my_filter_code2;
    x_filtered = filter(H, x1);

x_filtered为已经被滤波的信号，H为所设计的滤波器。
滤波效果如下：

由图可见，低频部分的分量被很好滤除，且对比该图和原小提琴信号的时域和频域图，时域部分看不出差别，频域部分除被滤除的部分，其他部分基本一致。最后将滤波之后的数据用audiowrite函数导出得到wav文件，其声音相对于混杂音频有极大的恢复且与原音频几乎没有差别。

附录：MATLAB代码

1、小提琴音频的加噪去噪处理代码

    clear;
clc;
format long;

% ==========原始信号========== %
[x1, fs1] = audioread('./实验二语音信号/小提琴.wav');

N1 = length(x1);  % 整个图由N1个样点构成
dt1 = 1 / fs1;
tscale1 = dt1 * N1;  % X轴显示的时间长度，单位为秒
t1 = 0 : dt1 : tscale1 - tscale1 / N1;

subplot(1, 2, 1);
plot(t1 .* 1000, x1);
title('小提琴信号时域图');
xlabel('t/ms', 'FontName', '宋体', 'FontWeight', 'normal', 'FontSize', 14);
ylabel('电压/V', 'FontName', '宋体', 'FontWeight', 'normal', 'FontSize', 14);
grid on;

y1 = fft(x1);
realy = 2 * abs(y1(1 : length(x1))) / length(x1);
realf = (0 : length(x1) - 1) * (fs1 / length(x1)); 
subplot(1, 2, 2);
stem(realf, realy, '.');
title('小提琴信号频谱图');
axis([0, 4000, 0, 0.04]);
xlabel('f/Hz', 'FontName', '宋体', 'FontWeight', 'normal', 'FontSize', 14);
ylabel('电压/V', 'FontName', '宋体', 'FontWeight', 'normal', 'FontSize', 14);
grid on;

figure;
% ==========噪声信号========== %
[x3, fs3] = audioread('./实验二语音信号/噪声.wav');

N3 = length(x3);  % 整个图由N1个样点构成
dt3 = 1 / fs3;
tscale3 = dt3 * N3;  % X轴显示的时间长度，单位为秒
t3 = 0 : dt3 : tscale3 - tscale3 / N3;

subplot(1, 2, 1);
plot(t3 .* 1000, x3);
title('噪声信号时域图');
xlabel('t/ms', 'FontName', '宋体', 'FontWeight', 'normal', 'FontSize', 14);
ylabel('电压/V', 'FontName', '宋体', 'FontWeight', 'normal', 'FontSize', 14);
grid on;

y3 = fft(x3);
realy = 2 * abs(y3(1 : length(x3))) / length(x3);
realf = (0 : length(x3) - 1) * (fs3 / length(x3)); 
subplot(1, 2, 2);
stem(realf, realy, '.');
title('噪声信号频谱图');
axis([0, 4000, 0, 0.04]);
xlabel('f/Hz', 'FontName', '宋体', 'FontWeight', 'normal', 'FontSize', 14);
ylabel('电压/V', 'FontName', '宋体', 'FontWeight', 'normal', 'FontSize', 14);
grid on;

figure;
% ==========小提琴混杂噪声信号========== %
[x2, fs2] = audioread('./实验二语音信号/小提琴混杂声音_缩混.wav');

N2 = length(x2);  % 整个图由N2个样点构成
dt2 = 1 / fs2;
tscale2 = dt2 * N2;  % X轴显示的时间长度，单位为秒
t2 = 0 : dt2 : tscale2 - tscale2 / N2;

subplot(1, 2, 1);
plot(t2 .* 1000, x2);
title('小提琴混杂噪声信号时域图');
xlabel('t/ms', 'FontName', '宋体', 'FontWeight', 'normal', 'FontSize', 14);
ylabel('电压/V', 'FontName', '宋体', 'FontWeight', 'normal', 'FontSize', 14);
grid on;

y2 = fft(x2);
realy = 2 * abs(y2(1 : length(x2))) / length(x2);
realf = (0 : length(x2) - 1) * (fs2 / length(x2)); 
subplot(1, 2, 2);
stem(realf, realy, '.');
title('小提琴混杂噪声信号频谱图');
axis([0, 4000, 0, 0.04]);
xlabel('f/Hz', 'FontName', '宋体', 'FontWeight', 'normal', 'FontSize', 14);
ylabel('电压/V', 'FontName', '宋体', 'FontWeight', 'normal', 'FontSize', 14);
grid on;

figure;
% ==========MATLAB工具箱生成的滤波器========== %
H = my_filter_code2;
x_filtered = filter(H, x1);

subplot(1, 2, 1);
plot(t1 .* 1000, x_filtered);
title('小提琴去噪信号时域图');
xlabel('t/ms', 'FontName', '宋体', 'FontWeight', 'normal', 'FontSize', 14);
ylabel('电压/V', 'FontName', '宋体', 'FontWeight', 'normal', 'FontSize', 14);
grid on;

y_filtered = fft(x_filtered);
realy = 2 * abs(y_filtered(1 : length(x_filtered))) / length(x_filtered);
realf = (0 : length(x_filtered) - 1) * (fs2 / length(x_filtered)); 
subplot(1, 2, 2);
stem(realf, realy, '.');
title('小提琴去噪信号频谱图');
axis([0, 4000, 0, 0.04]);
xlabel('f/Hz', 'FontName', '宋体', 'FontWeight', 'normal', 'FontSize', 14);
ylabel('电压/V', 'FontName', '宋体', 'FontWeight', 'normal', 'FontSize', 14);
grid on;

audiowrite('./实验二语音信号/小提琴去噪.wav', x_filtered, 8000);