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个滤波器的幅频、相频、衰减及每次滤波后的信号时域图。
![image-20211128145815600](https://i-blog.csdnimg.cn/blog_migrate/4b4cdab51bce21af78eadfeef15a46b2.png)
![image-20211128145835590](https://i-blog.csdnimg.cn/blog_migrate/27b9f8542739e9079934f77f868576e9.png)
![image-20211128145856791](https://i-blog.csdnimg.cn/blog_migrate/a40752ef58685043f5949f34bd543a56.png)
![image-20211128145920683](https://i-blog.csdnimg.cn/blog_migrate/6f5f24a1442888165a557bc2f5a17ad4.png)
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的最终输出信号的分析如下: