小波变换语音消噪(改进阈值)

改进阈值函数进行语音信号消噪，但是在程序运行过程中频频报错。本人经验不足调试不出，希望求得各位指导改进函数表达式附图

clear all; clc; close all;

fs=8000; %语音信号采样频率为8000

xx=wavread('lw1.wav');

x1=xx(:,1);%取单声道

t=(0:length(x1)-1)/8000;

y1=fft(x1,2048); %对信号做2048点FFT变换

f=fs*(0:1023)/2048;

figure(1)

plot(t,x1) %做原始语音信号的时域图形

y=awgn(x1',10,'measured'); %加10db的高斯白噪声

[snr,mse]=snrmse(x1,y')%求得信噪比均方误差

figure(2)

plot(t,y) %做加噪语音信号的时域图形

[c,l]=wavedec(y,3,'db1');%多尺度一维分解

%用db1小波对信号进行3层分解并提取系数

a3=appcoef(c,l,'db1',3);

%a2=appcoef(c,l,'db1',2);

%a1=appcoef(c,l,'db1',1);

d3=detcoef(c,l,3);

d2=detcoef(c,l,2);

d1=detcoef(c,l,1);

thr1=thselect(d1,'rigrsure');%阈值获取，使用Stein的无偏风险估计原理

thr2=thselect(d2,'rigrsure');

thr3=thselect(d3,'rigrsure');

%利用改进阈值函数进行去噪处理

gd1=Garrote_gg(d1,thr1);

gd2=Garrote_gg(d2,thr2);

gd3=Garrote_gg(d3,thr3);

c1=[a3 gd3 gd2 gd1];

y1=waverec(c2,l,'db1');%多尺度重构

[snr,mse]=snrmse(x1,y1')%求得信噪比均方误差

figure(3);

plot(t,y1);

function gd=Garrote_gg(a,b)%a为信号分解后的小波系数，b为获得的阈值

m=0.2*((a*a)-(b*b));

if (abs(a)>=b)

gd=sign(a)*(abs(a)-b/exp(m));

else (abs(a)

gd=0;

end

function [snr,mse]=snrmse(I,In)

% 计算信噪比函数

% I :原始信号

% In:去噪后信号

snr=0;

Ps=sum(sum((I-mean(mean(I))).^2));%signal power

Pn=sum(sum((I-In).^2)); %noise power

snr=10*log10(Ps/Pn);

mse=Pn/length(I);

QQ截图20130516175535.png(11.18 KB, 下载次数: 0)

改进函数表达式

本帖最后由罗志雄于2013-5-16 21:58 编辑

function [snr,mse]=snrmse(I,In)

% 计算信噪比函数

% I :原始信号

% In:去噪后信号

snr=0;

Ps=sum(sum((I-mean(mean(I))).^2));%signal power

Pn=sum(sum((I-In).^2)); %noise power

snr=10*log10(Ps/Pn);

mse=Pn/length(I);

修改后程序清单如下：

clear all; clc; close all;

fs=8000; %语音信号采样频率为8000

xx=wavread('lw1.wav');

x1=xx(:,1);%取单声道

x1=x1-mean(x1);

t=(0:length(x1)-1)/8000;

y1=fft(x1,2048); %对信号做2048点FFT变换

f=fs*(0:1023)/2048;

figure(1)

plot(t,x1) %做原始语音信号的时域图形

y=awgn(x1',10,'measured'); %加10db的高斯白噪声

[snr,mse]=snrmsel(x1',y) %求得信噪比均方误差

snr1=SNR_singlech(x1',y)

figure(2)

plot(t,y) %做加噪语音信号的时域图形

[c,l]=wavedec(y,3,'db1');%多尺度一维分解

%用db1小波对信号进行3层分解并提取系数

a3=appcoef(c,l,'db1',3);

%a2=appcoef(c,l,'db1',2);

%a1=appcoef(c,l,'db1',1);

d3=detcoef(c,l,3);

d2=detcoef(c,l,2);

d1=detcoef(c,l,1);

thr1=thselect(d1,'rigrsure');%阈值获取，使用Stein的无偏风险估计原理thr2=thselect(d2,'rigrsure');

thr3=thselect(d3,'rigrsure');

%利用改进阈值函数进行去噪处理

gd1=Garrote_gg(d1,thr1);

gd2=Garrote_gg(d2,thr2);

gd3=Garrote_gg(d3,thr3);

c1=[a3 gd3 gd2 gd1];

function gd=Garrote_gg(a,b)%a为信号分解后的小波系数，b为获得的阈值m=0.2*((a.*a)-(b*b));

if (abs(a)>=b)

gd=sign(a)*(abs(a)-b/exp(m));

else

gd=zeros(size(a));

end

y1=waverec(c1,l,'db1');%多尺度重构

[snr,mse]=snrmsel(x1',y1) %求得信噪比均方误差

figure(3);

plot(t,y1);