音频小波降噪实例（MATLAB）DWT-VoiceAnalysis

close all; clear; clc; % 实例功能： % 读取语音信号：'ReferAudio.flac' 并添加幅度为 NoiseAmplitude 的 gauss 白噪声。 % 采用 WaveName 小波对信号进行 level 尺度分解 % 采用 4 种方法对信号进行降噪处理并重构，计算相应的信噪比、均方根误差等对性能进行了比较。 % 保存的音频文件位于 AudioFile 内 % 最后修改时间：2022/11/18 %==============================================================author：Chen %% 参考信号定义 Num_plot = 1; % 图像编号 voice_sign = 0; % 音频播放标志 voice_sign = 0 时不显示是否播放音频。 NoiseAmplitude = 0.025;% 噪声强度 %========================================================================== % 音频文件名称 filename = 'ReferAudio.flac'; % 读取音频文件 [data_Audio, fs] = audioread(filename); % 关闭MATLAB警告 warning ('off') % 参数计算 St = (data_Audio(:,1))'; % 采用 data_Audio 1通道数据 N = length(St); % 信号长度 T = N/fs; % 信号持续时间 % 坐标轴初始化 t = linspace(0,T,N); f = linspace(-fs/2,fs/2,N); % 播放原声 if voice_sign if input('是否播放原声。 输入1播放\n') == 1 % 输入 1 播放原声 sound(St, fs) pause(T); % 等待音频播放结束 fprintf('播放完毕。\n') end end % 原信号频域 St_fft = 2*fft(St.*(hamming(N))')/N; % 对原信号×汉明窗来减小信号非周期截断产生的旁瓣 % 信号图像绘制 Fig_name = sprintf("信号时频图 Num %d",Num_plot); Num_plot = Num_plot+1; figure(NumberTitle="off",Name=Fig_name) subplot(2,1,1) plot(t,real(St)); xlabel('时间 s'); ylabel('幅度'); title('信号时域'); subplot(2,1,2) plot(f,fftshift(abs(St_fft))); xlabel('频率 Hz'); ylabel('幅度'); title('信号频域'); %% 添加噪声 gauss = NoiseAmplitude*randn(1,N); % 添加高斯白噪声 Yt = St + gauss; Yt_fft = 2*fft(Yt)/N; % 含噪声信号频域 % 保存含 NoiseAmplitude 强度的高斯白噪声 SaveFile_name = sprintf("./AudioFile/Gauss_%s.flac",extractBefore(filename,'.')); audiowrite(SaveFile_name,Yt,fs); % 含噪声信号绘制 Fig_name = sprintf("含噪声音频信号 Num %d",Num_plot); Num_plot = Num_plot+1; figure(NumberTitle="off",Name=Fig_name) subplot(2,1,1) % 含噪声参考信号时域 plot(t,real(Yt)); xlabel('时间 s'); ylabel('幅度'); title('含噪声音频信号时域'); subplot(2,1,2) % 含噪声参考信号频域 plot(f,fftshift(abs(Yt_fft))); xlabel('频率 Hz'); ylabel('幅度'); title('含噪声音频信号频域'); % 播放含燥声音频 if voice_sign if input('是否播放含噪声信号。 输入1播放\n') == 1 %1 % 输入 1 播放含燥声音频 sound(Yt, fs) pause(T); % 等待音频播放结束 fprintf('播放完毕。\n') end end %% 小波分解 % 参数定义================================================================== level = 3; WaveName = 'db5'; % 注：dbN （Daubechies小波 N 表示阶数）db1等同于Haar小波 % 多级分解时需要注意下采样是否满足奈奎斯特的问题； % Morlet 小波不能做离散小波变换和正交小波变换（不具备正交性，仅满足连续小波的允许条件） %========================================================================== [C, L] = wavedec(Yt, level, WaveName); % 获得尺度 level 的细节系数 for i = 1:level eval(['cd' num2str(i) '= detcoef(C, L, i)', ';']); % 命名为 cdx (x = 1:level) end % 获得尺度 level 的近似系数 eval(['ca' num2str(level) '= detcoef(C, L, level)', ';']); % 命名为 cax (x = level) % 小波分解图像绘制 Fig_name = sprintf("小波分解时域图像绘制 Num %d",Num_plot); Num_plot = Num_plot+1; figure(NumberTitle="off",Name=Fig_name) % 时域图像绘制 for i = 1:level % 细节系数时域绘制 subplot(level+1,1,i) plot(1:L(level +2-i),real(eval(['cd' num2str(i)]))); title(sprintf('细节系数：cd%d 时域',i)); end subplot(level+1,1,level+1) % 近似系数时域绘制 plot(1:L(1),real(eval(['ca' num2str(level)]))); title(sprintf('近似系数：ca%d 时域',level)); Fig_name = sprintf("小波分解频域图像绘制 Num %d",Num_plot); Num_plot = Num_plot+1; figure(NumberTitle="off",Name=Fig_name) % 频域图像绘制 % 细节系数频域绘制 for i = 1:level process_signal = eval(['cd' num2str(i)]); % 待处理信号 N_P = length(process_signal); % 待处理信号长度 % cDx 真实频域（x = 1:level） fft_ProcSign = 2*fft(process_signal)/N_P; subplot(level+1,1,i) %fw = linspace(-L(level +2-i)/2,L(level +2-i)/2,L(level +2-i)/fs*N); fw = (0:N_P/2)*(fs/N); % 绘制正频部分 plot(fw,abs(fft_ProcSign(1:N_P/2+1))); title(sprintf('cd%d 频域',i)); end % 近似系数频域 process_signal = eval(['ca' num2str(level)]); % 待处理信号 N_P = length(process_signal); % 待处理信号长度 % cAx 真实频域（x = level） fft_ProcSign = 2*fft(process_signal)/N_P; fw = (0:N_P/2)*(fs/N); % 绘制正频部分 subplot(level+1,1,level+1) plot(fw,abs(fft_ProcSign(1:N_P/2+1))); title(sprintf('近似系数：ca%d 频域',level)); %% 小波重构 （未降噪信号重构） s_rec = waverec(C, L, WaveName); Fig_name = sprintf("小波重构 Num %d",Num_plot); Num_plot = Num_plot+1; figure(NumberTitle="off",Name=Fig_name) subplot(2,1,1) plot(t,real(s_rec)) % 小波重构信号 xlabel('时间 s'); ylabel('幅度'); title('小波重构') subplot(2,1,2) plot(t,abs(s_rec - Yt)) % 重构信号与原信号差值 xlabel('时间 s'); ylabel('幅度'); title('重构信号与原信号差值') % 播放小波重构音频 if voice_sign if input('是否播放小波重构信号。 输入1播放\n') == 1 % 输入 1 播放小波重构信号 sound(s_rec, fs) pause(T); % 等待音频播放结束 fprintf('播放完毕。\n') end end %% 小波降噪 % 阈值模型采用软阈值降噪 method_num = 4; % 采用的方法数 method_name = ["缺省阈值", "penalty阈值", "细节系数百分比置零", "基于Stein分层阈值"];% 模型名称 %========================================================================== % 缺省的阈值模型 method 1 [thr1, sorh, keepapp] = ddencmp('den', 'wv', Yt); % 获得缺省阈值模型阈值 % thr1 = sqrt(2*log10(N)*sigma1); % 重建噪声信号 % 计算的阈值进行全局降噪（将th1应用与全部细节分量） xd1 = wdencmp('gbl', C, L, WaveName, level, thr1, 's', 1); xd1_fft = 2*fft(xd1)/N; % penalty阈值模型 method 2 sigma = wnoisest(C, L, 1); % 通过第1层细节系数估算信号的噪声强度 alpha = 2; % 选择参数：α = 2 thr2 = wbmpen(C, L, sigma, alpha); % 计算 Penalty阈值系数 xd2 = wdencmp('gbl', C, L, WaveName, level, thr2, 's', 1); xd2_fft = 2*fft(xd2)/N; % 通过将细节系数百分比设置为 0 来抑制噪声 method 3 Compress_DetLev = 1:level; % 初始化细节级别矩阵 Compress_Prop = zeros(1,level);% 初始化细节置零百分比矩阵 % 1:level-2 级别 90% 置零； level-1:level 级别 50% 置零。 Compress_Prop(1:level-2) = 99; Compress_Prop(level-1:level) = 80; xd3 = wthcoef("d",C,L,Compress_DetLev,Compress_Prop);% 得到新的矩阵C xd3 = waverec(xd3,L,WaveName); % 小波重构 xd3_fft = 2*fft(xd3)/N; % 采用分层思想：基于stein无偏估计的方法 method 4 [thr3, nkeep] = wdcbm(C, L, 1.001); % wdcbm(C,L,ALPHA) ALPHA 定义策略 [xd4, cxd, lxd, perf0, perfl2] = wdencmp('lvd', C, L, WaveName, level, thr3, 's'); xd4_fft = 2*fft(xd4)/N; %========================================================================== % method_num 种重构方法图像绘制 Fig_name = spri