Matlab语音信号分析与处理

%% 初始语音信号 [y0,fs0]=audioread('语音信号.wav'); %读取语音信号，fs为采样频率 n0=length(y0); figure('name','初始语音信号'); %绘出时域波形 x0=(0:1/n0:(1-1/n0))*10; subplot(2,1,1); plot(x0,y0); xlabel('时间/s'); ylabel('幅度'); title('时域波形'); %绘出频谱图 y0_fft=fft(y0,n0); %快速傅里叶变换 y0_Ayy=(abs(y0_fft)); %换算成实际的幅度 y0_Ayy=y0_Ayy/(n0/2); y0_Ayy(1)=y0_Ayy(1)/2; y0_fn=1:n0; for i=1:n0 %换算成实际的频率值 y0_fn(i) = (i-1)*fs0/n0; end subplot(2,1,2); plot(y0_fn(1:n0/2),y0_Ayy(1:n0/2)); xlabel('频率/Hz'); ylabel('幅度'); title('频谱图'); % %绘出相频特性曲线 % Pyy=1:n/2; % for i=1:n/2 % Pyy(i)=atan2(imag(y1(i)),real(y1(i))); %计算相位 % Pyy(i)=Pyy(i)*180/pi; %换算为角度 % end % subplot(2,2,4); % plot(fn(1:n/2),Pyy(1:n/2)); %显示相位图 % title('相频特性曲线'); %% 信号降采样处理 fs=21000; y1=resample(y0,fs,fs0); %采样从fs0降到fs n1=length(y1); figure('name','采样后的语音信号'); %绘出时域波形 x1=(0:1/n1:(1-1/n1))*10; subplot(2,1,1); plot(x1,y1); xlabel('时间/s'); ylabel('幅度'); title('时域波形'); %绘出频谱图 y1_fft=fft(y1,n1); %快速傅里叶变换 y1_Ayy=(abs(y1_fft)); %换算成实际的幅度 y1_Ayy=y1_Ayy/(n1/2); y1_Ayy(1)=y1_Ayy(1)/2; y1_fn=1:n1; for i=1:n1 %换算成实际的频率值 y1_fn(i) = (i-1)*fs/n1; end subplot(2,1,2); plot(y1_fn(1:n1/2),y1_Ayy(1:n1/2)); xlabel('频率/Hz'); ylabel('幅度'); title('频谱图'); %% 加高斯白噪声 y2=awgn(y1,30); n2=length(y2); figure('name','加噪声后的语音信号'); %绘出时域波形 x2=(0:1/n2:(1-1/n2))*10; subplot(2,1,1); plot(x2,y2); xlabel('时间/s'); ylabel('幅度'); title('时域波形'); %绘出频谱图 y2_fft=fft(y2,n2); %快速傅里叶变换 y2_Ayy=(abs(y2_fft)); %换算成实际的幅度 y2_Ayy=y2_Ayy/(n2/2); y2_Ayy(1)=y2_Ayy(1)/2; y2_fn=1:n2; for i=1:n2 %换算成实际的频率值 y2_fn(i) = (i-1)*fs/n2; end subplot(2,1,2); plot(y2_fn(1:n2/2),y2_Ayy(1:n2/2)); xlabel('频率/Hz'); ylabel('幅度'); title('频谱图'); %% 对加噪声后的信号进行调制 t=(0:1/n2:(1-1/n2))*10; y3=(1+y2).*cos(12000*pi*t'); n3=length(y3); figure('name','振幅调制(AM)'); %绘出时域波形 subplot(3,1,1); plot(t,y3); xlabel('时间/s'); ylabel('幅度'); title('时域波形'); %绘出频谱图 y3_fft=fft(y3,n3); %快速傅里叶变换 y3_Ayy=(abs(y3_fft)); %换算成实际的幅度 y3_Ayy=y3_Ayy/(n3/2); y3_Ayy(1)=y3_Ayy(1)/2; y3_fn=1:n3; for i=1:n3 %换算成实际的频率值 y3_fn(i) = (i-1)*fs/n3; end subplot(3,1,2); plot(y3_fn(1:n3/2),y3_Ayy(1:n3/2)); xlabel('频率/Hz'); ylabel('幅度'); title('频谱图1'); subplot(3,1,3); plot(y3_fn(1:n3/2),y3_Ayy(1:n3/2)); axis([0,12000,0,0.004]); xlabel('频率/Hz'); ylabel('幅度'); title('频谱图2'); %% 理想带通传输信道传输 y4=bandpass(y3,4000,8000,fs); n4=length(y4); figure('name','通过理想带通传输信道的信号'); %绘出时域波形 x4=(0:1/n4:(1-1/n4))*10; subplot(3,1,1); plot(x4,y4); xlabel('时间/s'); ylabel('幅度'); title('时域波形'); %绘出频谱图 y4_fft=fft(y4,n4); %快速傅里叶变换 y4_Ayy=(abs(y4_fft)); %换算成实际的幅度 y4_Ayy=y4_Ayy/(n4/2); y4_Ayy(1)=y4_Ayy(1)/2; y4_fn=1:n4; for i=1:n4 %换算成实际的频率值 y4_fn(i) = (i-1)*fs/n4; end subplot(3,1,2); plot(y4_fn(1:n4/2),y4_Ayy(1:n4/2)); xlabel('频率/Hz'); ylabel('幅度'); title('频谱图'); subplot(3,1,3); plot(y4_fn(1:n4/2),y4_Ayy(1:n4/2)); axis([0,12000,0,0.004]); xlabel('频率/Hz'); ylabel('幅度'); title('频谱图'); %% 对接收到的信号进行解调 y5=envelope(y4)-1; n5=length(y5); figure('name','解调后的信号'); %绘出时域波形 x5=(0:1/n5:(1-1/n5))*10; subplot(3,1,1); plot(x5,y5); axis([0,10,-0.6,0.4]); xlabel('时间/s'); ylabel('幅度'); title('时域波形'); %绘出频谱图 y5_fft=fft(y5,n5); %快速傅里叶变换 y5_Ayy=(abs(y5_fft)); %换算成实际的幅度 y5_Ayy=y5_Ayy/(n5/2); y5_Ayy(1)=y5_Ayy(1)/2; y5_fn=1:n5; for i=1:n5 %换算成实际的频率值 y5_fn(i) = (i-1)*fs/n5; end subplot(3,1,2); plot(y5_fn(1:n5/2),y5_Ayy(1:n5/2)); xlabel('频率/Hz'); ylabel('幅度'); title('频谱图1'); subplot(3,1,3); plot(y5_fn(1:n5/2),y5_Ayy(1:n5/2)); axis([0,12000,0,0.004]); xlabel('频率/Hz'); ylabel('幅度'); title('频谱图2'); %% 去噪 y6=wiener(y5,y1); y6=wiener(y6,y1); y6=filter(lowpass,y6); n6=length(y6); figure('name','去噪后的信号'); %绘出时域波形 x6=(0:1/n6:(1-1/n6))*10; subplot(2,1,1); plot(x6,y6); xlabel('时间/s'); ylabel('幅度'); title('时域波形'); %绘出频谱图 y6_fft=fft(y6,n6); %快速傅里叶变换 y6_Ayy=(abs(y6_fft)); %换算成实际的幅度 y6_Ayy=y6_Ayy/(n6/2); y6_Ayy(1)=y6_Ayy(1)/2; y6_fn=1:n6; for i=1:n6 %换算成实际的频率值 y6_fn(i) = (i-1)*fs/n6; end subplot(2,1,2); plot(y6_fn(1:n6/2),y6_Ayy(1:n6/2)); xlabel('频率/Hz'); ylabel('幅度'); title('频谱图');