MATLAB读取声音信号进行频域阈值处理

clc
[pyr,fs]=audioread('guo.wav.m4a');%声音读取
n=length(pyr);
pyr1=fft(pyr,n);                            %快速傅里叶变换
figure(1)
subplot(2,1,1);
plot(pyr);   
xlabel('时间');
ylabel('幅度');
title('初始信号波形');   %绘出时域波
grid
subplot(2,1,2);                         %绘出频域频谱
plot(abs(fftshift(pyr1)));
title('初始信号频谱');
xlabel('频率');
ylabel('幅度');
grid
noise=0.1*randn(n,1);%加噪声
s=pyr+noise;
pause;
n=length(s);   %画出加噪之后，其时域频域
S=fft(s,n);
figure(2)
subplot(2,1,1)
plot(s);
title('加噪声后信号波形')
xlabel('时间');
ylabel('幅度');
grid;
subplot(2,1,2)
plot(abs(fftshift(S)));
xlabel('频率');
ylabel('幅度');
title('加噪声后信号信号频谱');
grid;
pause
[c,l]=wavedec(n,3,'db1');
 a3=appcoef(c,l,'db1',3);
 d3=detcoef(c,l,3);
 d2=detcoef(c,l,2);
 d1=detcoef(c,l,1);
 thr=1;%硬阈值处理
 ythard1=wthresh(d1,'h',thr);
 ythard2=wthresh(d2,'h',thr);
 ythard3=wthresh(d3,'h',thr);
 c2=[a3 ythard3 ythard2 ythard1];
 s3=waverec(c2,l,'db1');
  %软阈值处理
 ytsofrd1=wthresh(d1,'s',thr);
 ytsofrd2=wthresh(d2,'s',thr);
 ytsofrd3=wthresh(d3,'s',thr);
 c3=[a3 ytsofrd3 ytsofrd2 ytsofrd1];
 s4=waverec(c3,l,'db1'); 
 subplot(2,1,1);plot(s3);title('硬阈值处理');
 subplot(2,1,2);plot(s4);title('软阈值处理');