matlab小波去噪

Question

matlab进行小波去噪的步骤是哪些？第一步，信号分解，采用wavedec函数。设置分解数N和小波基函数。第二步，阈值处理，采用软
Donoho-Johnstone阈值应用于小波系数，并使用电平噪声的电平估计进行重新缩放。第二步怎么实现呢？

Answer 1

wavedec小波分解，[thrl,sorh,keepapp]=ddencmp('den','wv',y)， s2=wdencmp('gbl',c1,l,'db4',4,thrl,sorh,keepapp)软阈值处理：

clc,clear
[x,Fs]= audioread('MUsic_Test.wav');
snr = 10;  %设定信噪比，单位db
noise = randn(size(x));         % 用randn函数产生高斯白噪声
Nx = length(x);                 % 求出信号x长
signal_power = 1/Nx*sum(x.*x);  % 求出信号的平均能量
noise_power = 1/Nx*sum(noise.*noise);  % 求出噪声的能量
noise_variance = signal_power / ( 10^(snr/10) );  % 计算出噪声设定的方差值
noise = sqrt(noise_variance/noise_power)*noise;   % 按噪声的平均能量构成相应的白噪声
y = x + noise;                  % 合成带噪语音
[c,l]=wavedec(y,10,'db4');  %db4小波3级分解  https://ww2.mathworks.cn/help/wavelet/ref/wavedec.html?s_tid=doc_ta
ca3=appcoef(c,l,'db4',3)';
cd3=detcoef(c,l,3);
cd2=detcoef(c,l,2);
cd1=detcoef(c,l,1);

cdd3=zeros(1,length(cd3));
cdd2=zeros(1,length(cd2));
cdd1=zeros(1,length(cd1));
c1=[ca3,cdd3,cdd2,cdd1];
[thrl,sorh,keepapp]=ddencmp('den','wv',y);  %返回小波阈值去噪的默认参数-[thr,sorh,keepapp] = ddencmp(in1,in2,x) ...
% returns default values for denoising or compression, using wavelets or wavelet packets, of the input data x. x is a real-valued vector ...
% or 2-D matrix. thr is the threshold, and sorh indicates soft or hard thresholding. ...
% keepapp can be used as a flag to set whether or not the approximation coefficients are thresholded.
s2=wdencmp('gbl',c1,l,'db4',4,thrl,sorh,keepapp);  %小波重构
%绘图
windowLength = 256;  %帧长，语谱图参数
win = hamming(windowLength,'periodic');  %窗口函数（汉明窗）
overlap = 128; %帧移（一般为帧长的一半）
ffTLength = windowLength; %做DFT的点数，一般和帧长一样
subplot(3,2,1)
plot(x)
subplot(3,2,2)
spectrogram(x,win,overlap,ffTLength,Fs,'yaxis')
subplot(3,2,3) 
plot(y)
subplot(3,2,4)
spectrogram(y,win,overlap,ffTLength,Fs,'yaxis'); 
subplot(3,2,5)
plot(s2);
subplot(3,2,6)
spectrogram(s2,win,overlap,ffTLength,Fs,'yaxis'); 
sound(s2,Fs)
signal_power = sum(abs(x').^2)/length(x');
noise_Power0 = sum(abs(y'-x').^2)/length(y'-x');
SNR0 = 10*log10(signal_power/noise_Power0);
disp(['带噪音频信噪比：',num2str(SNR0)])
noise_Power2 = sum(abs(s2-x').^2)/length(s2-x');  %（带噪信号-纯信号）的平方
SNR = 10*log10(signal_power/noise_Power2);
disp(['强制消噪后信噪比：',num2str(SNR)])