Matlab 音频信号处理

Posted 2022-06-15 最忆是江南.

Matlab 音频信号处理

  现有原始音频文件SunshineSquare.wav，后半段音频被人为加上了多频段的噪音，使用matlab工具对其进行分析并消除噪音还原出无噪音频文件。

 

Step1：分析原始音频文件

读取

[audio_data, fs] = audioread('../../SunshineSquare.wav');
L = length(audio_data);

% 听一听该音频
% soundsc(audio_data, fs);

绘制时域图

  从下面时域图中可以看出，源音频文件时长11秒左右，从第8秒左右开始为一段强噪声。

figure('Name','Original');


tt = (1 : L) * 1/fs;
subplot(4, 1, 1);plot(tt, audio_data);grid;
title('Time Domain');xlabel('时间（s）');ylabel('Amp');

频谱分析

  使用快速傅里叶变换FFT对时域信号进行频谱分析，绘制的频域图反映出有四段噪声频段，后续消除噪声工作就是对这四个频段进行滤波。

L = 2^nextpow2(L); % 先从原始信号长度确定下一个 2 次幂的新长度,用尾随零填充信号以改善 fft 的性能
Spectrum = fft(audio_data,L); % 快速傅里叶变换，转换为频域
f_x = fs*(0:(L/2))/L; % 定义频域，即变换横坐标量纲为频率
mag_FFT_audio_data = abs(Spectrum / L); % 取幅值，并归一化
ang_FFT_audio_data = angle(Spectrum)*180/pi;

subplot(4, 1, 2);plot(f_x,mag_FFT_audio_data(1 : L/2+1)) 
title('Frequency Domain');xlabel('频率（Hz）');ylabel('mag');
subplot(4, 1, 3);plot(f_x,ang_FFT_audio_data(1 : L/2+1)) 
title('Frequency Domain');xlabel('频率（Hz）');ylabel('ang');

绘制声谱图

  声谱图能够以颜色反映各频段信号强度随时间的变化和分布情况，如图所示，从第8秒左右开始有4个频段的强度持续较高，这就是上面频域分析出的4个噪声频段。但需要注意到，在噪声的开始和结束时刻，所有频段的强度都很高，在图中显示为两条垂直的黄线，这两条线仅靠滤波无法消除。

subplot(4, 1, 4)
spectrogram(audio_data,hann(256),250,256,fs,'yaxis');

 

 

Step2：滤波

滤波法1——手动衰减

  手动衰减各个噪声频段的幅度，注意FFT变换后的频谱图为偶对称，所以共有4对（8个）频段，逐一消除。

[m, n] = max(Spectrum)
Spectrum(n-850: n+850) = 0; % 第一个噪声频段

[m, n] = max(Spectrum)
Spectrum(n-850: n+850) = 0; % 第二个噪声频段

[m, n] = max(Spectrum)
Spectrum(n-850: n+850) = 0; % 第三个噪声频段

[m, n] = max(Spectrum)
Spectrum(n-850: n+850) = 0; % 第四个噪声频段

[m, n] = max(Spectrum)
Spectrum(n-850: n+850) = 0; % 第五个噪声频段

[m, n] = max(Spectrum)
Spectrum(n-850: n+850) = 0; % 第六个噪声频段

[m, n] = max(Spectrum)
Spectrum(n: n+850) = 0; % 第七个噪声频段

[m, n] = max(Spectrum)
Spectrum(n-850: n+850) = 0; % 第八个噪声频段

% 逆fft获得处理后的音频
ifft_data = real(ifft(Spectrum, L));
audio_data = ifft_data(1:length(audio_data));

滤波效果

  绘制滤波后音频信号的时域图、频域图、声谱图，可看见时域图上噪声被消除，频谱图上4个频段的信号全部被消除，声谱图有四条水平的消磨痕迹。

 

滤波法2——数字滤波器

  设计梳状滤波器，因为有4个噪声频段，所以设计4个带阻滤波器逐一滤波。

%h1[n] = [1,-2,1];
h1 = [1,-2,1];
h2 = [1,-1.24698,1];
h3 = [1,0.44504,1];
h4 = [1,1.80194,1];
a = [1,0,0];

%滤波过程
audio_data_h1 = filter(h1,a,audio_data);
audio_data_h2 = filter(h2,a,audio_data_h1);
audio_data_h3 = filter(h3,a,audio_data_h2);
audio_data_h4 = filter(h4,a,audio_data_h3);

audio_data = audio_data_h4;

滤波器分析

  绘制4个滤波器的幅频、相频、衰减及每次滤波后的信号时域图。

 

 

Step3：后续处理

  上面已分析，在噪声开始和结束时刻都有全频段的高强度信号，不是滤波能消除的，这里采用时域衰减，即在时域采用同第一种滤波方法相同的思想，手动衰减这两个时刻的信号强度。随后输出处理过的音频文件。

[m, n] = max(audio_data);
audio_data(n-200: n+200) = 0.0001 * audio_data(n-200: n+200); % 衰减第一个高音时段
[m, n] = max(audio_data);
audio_data(n-200: n+200) = 0.0001 * audio_data(n-200: n+200); % 衰减第二个高音频段

% 听一下处理后的音频
% soundsc(audio_data, fs);

% 输出音频文件
audiowrite('SunshineSquare_Processed.wav',audio_data,fs); 

使用滤波法2的最终输出信号的分析如下：

开发者涨薪指南 48位大咖的思考法则、工作方式、逻辑体系