用MATLAB对语音信号进行小波分解重构

suffer

用MATLAB对一语音信号进行小波分解，然后对其各层系数进行处理以达到小波抑制的目的，重构处理后的信号，画出波形分析。

1. %装载原始信号
   
2. load sumsin;
   
3. s=sumsin;
   
4. %==============================
   
5. %设置小波名并利用coif3小波进行4层分解
   
6. w='coif3';
   
7. maxlev=4;
   
8. [c,l]=wavedec(s,maxlev,w);
   
9. newc=c;
   
10. %==============================
    
11. %将分解后的第三、四层细节系数值为0
    
12. newc=wthcoef('d',c,l,[3,4]);
    
13. %==============================
    
14. %在原始信号的时间区间[400，600]内将第一层细节系数值为0
    
15. %并且将其他系数进行衰减，求出第一层系数起始点和终止点的
    
16. %索引值
    
17. k=maxlev+1;
    
18. first=sum(l(1:k-1))+1;
    
19. last=first+l(k-1);
    
20. indd1=first:last;
    
21. %==============================
    
22. %将系数除以3，进行信号衰减
    
23. newc(indd1)=c(indd1)/3;
    
24. %==============================
    
25. %在区间[400，600]上求出第一层系数索引
    
26. indd1=(first+400/2):(first+600/2);
    
27. %==============================
    
28. %将该索引值置为0
    
29. newc(indd1)=zeros(size(indd1));
    
30. %==============================
    
31. %将第二层中相应于原始信号t=500的时间点处的系数置为4
    
32. k=maxlev;
    
33. first=sum(l(1:k-1))+1;
    
34. newc(first+500/2^2)=4;
    
35. %==============================
    
36. %综合修改后的分解结构
    
37. synth=waverec(newc,l,w);
    
38. %==============================
    
39. %用图示出上述修改结果
    
40. subplot(2,2,1);
    
41. plot(s);
    
42. title('原始信号');
    
43. subplot(2,2,2);
    
44. plot(c);
    
45. title('coif3小波分解后的系数');
    
46. subplot(2,2,3);
    
47. plot(synth);
    
48. title('小波抑制后的信号');
    
49. subplot(2,2,4);
    
50. plot(newc);
    
51. title('修改后的小波分解系数');

复制代码

suffer

用MATLAB实现对一语音信号用不同小波进行分解，提取各层的高低频系数，画出各系数波形，并重构。

1. t=0:1:100*pi;
   
2. s=sin(3*t)+sin(0.3t)+sin(0.03t);
   
3. subplot(6,2,1);plot(s);
   
4. title('原始信号s');
   
5. %====================================
   
6. %对s进行小波分解：db3   5层
   
7. [c,l]=wavedec(s,5,'db3');
   
8. %====================================
   
9. %提取小波分解的低频系数
   
10. a5=appcoef(c,l,'db3',5);
    
11. a4=appcoef(c,l,'db3',4);
    
12. a3=appcoef(c,l,'db3',3);
    
13. a2=appcoef(c,l,'db3',2);
    
14. a1=appcoef(c,l,'db3',1);
    
15. %====================================
    
16. %提取小波分解的各层高频系数
    
17. d5=detcoef(c,l,5);
    
18. d4=detcoef(c,l,4);
    
19. d3=detcoef(c,l,3);
    
20. d2=detcoef(c,l,2);
    
21. d1=detcoef(c,l,1);       
    
22. %====================================
    
23. %绘出各系数的图形
    
24. subplot(6,2,3);plot(a5);
    
25. Ylabel('a5');
    
26. subplot(6,2,5);plot(a4);
    
27. Ylabel('a4');
    
28. subplot(6,2,7);plot(a3);
    
29. Ylabel('a3');
    
30. subplot(6,2,9);plot(a2);
    
31. Ylabel('a2');
    
32. subplot(6,2,11);plot(a1);
    
33. Ylabel('a1');
    
34. subplot(6,2,4);plot(d5);
    
35. Ylabel('d5');
    
36. subplot(6,2,6);plot(d4);
    
37. Ylabel('d4');
    
38. subplot(6,2,6);plot(d3);
    
39. Ylabel('d3');
    
40. subplot(6,2,8);plot(d2);
    
41. Ylabel('d2');
    
42. %====================================
    
43. %重构信号s
    
44. s1=waverec(c,l,'db1');
    
45. subplot(5,2,9);plot(s1);
    
46. Ylabel('s1');
    
47. %====================================
    
48. %下面用小波‘coif3’重复上述过程
    
49. [c,l]=wavedec(s,3,'coif3');
    
50. a3=appcoef(c,l,'coif3',3);
    
51. d3=detcoef(c,l,3);
    
52. d2=detcoef(c,l,2);
    
53. d1=detcoef(c,l,1);
    
54. subplot(5,2,2);plot(a3);
    
55. Ylabel('a3');
    
56. subplot(5,2,4);plot(d3);
    
57. Ylabel('d3');
    
58. subplot(5,2,6);plot(d2);
    
59. Ylabel('d2');
    
60. subplot(5,2,8);plot(d1);
    
61. Ylabel('d1');
    
62. s2=waverec(c,l,'coif3');
    
63. subplot(5,2,10);plot(s2);
    
64. Ylabel('s2');

复制代码

suffer

用MATLAB对一图象分别用不同小波分解，观察高低频部分进行重构，比较重构误差，再进行阈值消噪，比较消噪前后图象。

1. load noiswom; %装载图像信号
   
2. whos
   
3. [swa,swh,swv,swd] = swt2(X,1,'db1'); %完成图像的单层次小波分解
   
4. whos %观察各个系数向量的结构
   
5. figure(1);
   
6. map = pink(size(map,1)); %显示低频和高频系数
   
7. colormap(map)
   
8. subplot(2,2,1), image(wcodemat(swa,192));
   
9. title('Approximation swa')
   
10. subplot(2,2,2), image(wcodemat(swh,192));
    
11. title('Horiz. Detail swh')
    
12. subplot(2,2,3), image(wcodemat(swv,192));
    
13. title('Vertical Detail swv')
    
14. subplot(2,2,4), image(wcodemat(swd,192));
    
15. title('Diag. Detail swd')
    
16. A0 = iswt2(swa,swh,swv,swd,'db1'); %通过平稳小波逆变换重构图像
    
17. err = max(max(abs(X-A0))) %检查重构误差
    
18. nulcfs = zeros(size(swa)); %从第三步中产生的系数swa、swh、swv和swd构造第一层的低频和高频（A1、H1、V1和D1）部分
    
19. A1 = iswt2(swa,nulcfs,nulcfs,nulcfs,'db1');
    
20. H1 = iswt2(nulcfs,swh,nulcfs,nulcfs,'db1');
    
21. V1 = iswt2(nulcfs,nulcfs,swv,nulcfs,'db1');
    
22. D1 = iswt2(nulcfs,nulcfs,nulcfs,swd,'db1');
    
23. figure(2); %显示第一层的分解结果
    
24. colormap(map)
    
25. subplot(2,2,1), image(wcodemat(A1,192));
    
26. title('Approximation A1')
    
27. subplot(2,2,2), image(wcodemat(H1,192));
    
28. title('Horiz. Detail H1')
    
29. subplot(2,2,3), image(wcodemat(V1,192));
    
30. title('Vertical Detail V1')
    
31. subplot(2,2,4), image(wcodemat(D1,192));
    
32. title('Diag. Detail D1')
    
33. [swa,swh,swv,swd] = swt2(X,3,'db1'); %采用db1来完成多尺度二维离散平稳小波分解
    
34. clear A0 A1 D1 H1 V1 err nulcfs %观察swa,swh,swv,swd的存储结构
    
35. whos
    
36. figure(3); %显示多尺度二维平稳小波分解结果
    
37. colormap(map)
    
38. kp = 0;
    
39. for i = 1:3
    
40. subplot(3,4,kp+1), image(wcodemat(swa(:,:,i),192));
    
41. title(['Approx. cfs level ',num2str(i)])
    
42. subplot(3,4,kp+2), image(wcodemat(swh(:,:,i),192));
    
43. title(['Horiz. Det. cfs level ',num2str(i)])
    
44. subplot(3,4,kp+3), image(wcodemat(swv(:,:,i),192));
    
45. title(['Vert. Det. cfs level ',num2str(i)])
    
46. subplot(3,4,kp+4), image(wcodemat(swd(:,:,i),192));
    
47. title(['Diag. Det. cfs level ',num2str(i)])
    
48. kp = kp + 4;
    
49. end
    
50. mzero = zeros(size(swd)); %从系数中重构第3层的低频信号
    
51. A = mzero;
    
52. A(:,:,3) = iswt2(swa,mzero,mzero,mzero,'db1');
    
53. H = mzero; V = mzero; %重构第1、2、3层的高频信号
    
54. D = mzero;
    
55. for i = 1:3
    
56. swcfs = mzero; swcfs(:,:,i) = swh(:,:,i);
    
57. H(:,:,i) = iswt2(mzero,swcfs,mzero,mzero,'db1');
    
58. swcfs = mzero; swcfs(:,:,i) = swv(:,:,i);
    
59. V(:,:,i) = iswt2(mzero,mzero,swcfs,mzero,'db1');
    
60. swcfs = mzero; swcfs(:,:,i) = swd(:,:,i);
    
61. D(:,:,i) = iswt2(mzero,mzero,mzero,swcfs,'db1');
    
62. end
    
63. A(:,:,2) = A(:,:,3) + H(:,:,3) + V(:,:,3) + D(:,:,3); %通过第3层的低频部分和第1、2、3层的高频部分重构第1、2层的低频部分
    
64. A(:,:,1) = A(:,:,2) + H(:,:,2) + V(:,:,2) + D(:,:,2);
    
65. figure(4); %显示第1、2、3层的低频和高频部分
    
66. colormap(map)
    
67. kp = 0;
    
68. for i = 1:3
    
69. subplot(3,4,kp+1), image(wcodemat(A(:,:,i),192));
    
70. title(['Approx. level ',num2str(i)])
    
71. subplot(3,4,kp+2), image(wcodemat(H(:,:,i),192));
    
72. title(['Horiz. Det. level ',num2str(i)])
    
73. subplot(3,4,kp+3), image(wcodemat(V(:,:,i),192));
    
74. title(['Vert. Det. level ',num2str(i)])
    
75. subplot(3,4,kp+4), image(wcodemat(D(:,:,i),192));
    
76. title(['Diag. Det. level ',num2str(i)])
    
77. kp = kp + 4;
    
78. end;
    
79. thr = 44.5; %阈值消噪
    
80. sorh = 's';
    
81. dswh = wthresh(swh,sorh,thr);
    
82. dswv = wthresh(swv,sorh,thr);
    
83. dswd = wthresh(swd,sorh,thr);
    
84. clean = iswt2(swa,dswh,dswv,dswd,'db1');
    
85. thr = 44.5;
    
86. sorh = 's';
    
87. dswh = wthresh(swh,sorh,thr);
    
88. dswv = wthresh(swv,sorh,thr);
    
89. dswd = wthresh(swd,sorh,thr);
    
90. clean = iswt2(swa,dswh,dswv,dswd,'db1');
    
91. figure(5); %对比消噪前后的图像
    
92. colormap(map)
    
93. subplot(1,2,1), image(wcodemat(X,192));
    
94. title('Original image')
    
95. subplot(1,2,2), image(wcodemat(clean,192));
    
96. title('De-noised image')
    
97. end

复制代码

suffer

用MATLAB对一语音信号进行小波分解，分别用强阈值，软阈值，默认阈植进行消噪处理。

1. %装载采集的信号leleccum.mat
   
2. load leleccum;
   
3. %=============================
   
4. %将信号中第2000到第3450个采样点赋给s
   
5. indx=2000:3450;
   
6. s=leleccum(indx);
   
7. %=============================
   
8. %画出原始信号
   
9. subplot(2,2,1);
   
10. plot(s);
    
11. title('原始信号');
    
12. %=============================
    
13. %用db1小波对原始信号进行3层分解并提取系数
    
14. [c,l]=wavedec(s,3,'db1');
    
15. a3=appcoef(c,l,'db1',3);
    
16. d3=detcoef(c,l,3);
    
17. d2=detcoef(c,l,2);
    
18. d1=detcoef(c,l,1);
    
19. %=============================
    
20. %对信号进行强制性消噪处理并图示结果
    
21. dd3=zeros(1,length(d3));
    
22. dd2=zeros(1,length(d2));
    
23. dd1=zeros(1,length(d1));
    
24. c1=[a3 dd3 dd2 dd1];
    
25. s1=waverec(c1,l,'db1');
    
26. subplot(2,2,2);
    
27. plot(s1);grid;
    
28. title('强制消噪后的信号');
    
29. %=============================
    
30. %用默认阈值对信号进行消噪处理并图示结果
    
31.   %用ddencmp函数获得信号的默认阈值
    
32. [thr,sorh,keepapp]=ddencmp('den','wv',s);
    
33. s2=wdencmp('gbl',c,l,'db1',3,thr,sorh,keepapp);
    
34. subplot(2,2,3);
    
35. plot(s2);grid;
    
36. title('默认阈值消噪后的信号');
    
37. %=============================
    
38. %用给定的软阈值进行消噪处理
    
39. softd1=wthresh(d1,'s',1.465);
    
40. softd2=wthresh(d2,'s',1.823);
    
41. softd3=wthresh(d3,'s',2.768);
    
42. c2=[a3 softd3 softd2 softd1];
    
43. s3=waverec(c2,l,'db1');
    
44. subplot(2,2,4);
    
45. plot(s3);grid;
    
46. title('给定软阈值消噪后的信号');

复制代码

xiao1982

回复 suffer 的帖子

请教：softd1=wthresh(d1,'s',1.465);softd2=wthresh(d2,'s',1.823);softd3=wthresh(d3,'s',2.768);中1.465、1.823、2.768这三个数值是怎么确定的？谢谢！

jrjian

对c1=[a3 dd3 dd2 dd1];
不是很懂呢，希望解释下，谢谢

happy

xiao1982 发表于 2010-9-29 01:51

                               
登录/注册后可看大图

回复 suffer 的帖子

请教：softd1=wthresh(d1,'s',1.465);softd2=wthresh(d2,'s',1.823);softd3=wthresh(d ...

根据基本的噪声类型，阈值的选取有以下4个规则，其中每一条规则对应于函数thselect中输入参数tptr的一个选项。
（1）选项tptr='rigrsure'，是一种基于Stein无偏似然估计原理的自适应阈值选择。给定一个阈值t，得到它的似然估计，再将其似然最小化，就可得到所选的阈值。这是一种软件阈值估计器。
（2）选项tptr='sqtwolog'，是一种固定的阈值形式，它所产生的阈值为sqrt(2*log(length(X)))。
（3）选项tptr='heursure'，是前两种阈值的综合，所选择的是最优预测变量阈值。如果信噪比很小，而SURE估计有很大的噪声，此时就需要采用这种固定的阈值形式。
（4）选项tptr='minimaxi'，也是一种固定的阈值形式，它所产生的是一个最小均方差的极值，而不是无误差。

happy

jrjian 发表于 2010-10-24 23:37

                               
登录/注册后可看大图

对c1=[a3 dd3 dd2 dd1];
不是很懂呢，希望解释下，谢谢

这是强制消噪处理
强制消噪处理是把小波分解结构中的高频系数全部变为0,即把高频部分全部滤除掉,然后再对信号进行重构处理

jrjian

当处理自己录制的语音时，总是在c1=[a3 dd3 dd2 dd1];
地方出错，这几个参数都是小波分解自动产生的吧，出错的原因会在哪呢。谢谢

happy

jrjian 发表于 2010-10-26 11:08

                               
登录/注册后可看大图

当处理自己录制的语音时，总是在c1=[a3 dd3 dd2 dd1];
地方出错，这几个参数都是小波分解自动产生的吧，出 ...

估计会是数据长度等问题造成的
你也不给个具体的错误信息之类的，别人没办法帮你

另外你所说“这几个参数”指得是什么？

ximeiji

我在编写图像去噪时也总是在c1=[a3 dd3 dd2 dd1];
这里出现错误，好像是维数不对，请问要怎么改啊？很急，希望您能帮我看看！