爱华网基于MATLAB的语音信号处理
本文转载自龙行天下《基于MATLAB的语音信号处理》
程序1:用MATLAB对原始语音信号进行分析，画出它的时域波形和频谱fs=8000;%语音信号采样频率为8000
x1=wavread('pb8k.wav');
t=(0:length(x1)-1)/8000;
y1=fft(x1,2048);%对信号做2048点FFT变换
f=fs*(0:1023)/2048;
figure(1)
plot(t,x1)%做原始语音信号的时域图形
grid on;axis tight;
title('原始语音信号');
xlabel('time(s)');
ylabel('幅度');
figure(2)
plot(f,abs(y1(1:1024)))%做原始语音信号的FFT频谱图
grid on;axis tight;
title('原始语音信号FFT频谱')
xlabel('Hz');
ylabel('幅度');
**************************************************************************************************************************
程序2:给原始的语音信号加上一个高频余弦噪声，频率为3.8kHz。画出加噪后的语音信号时域和频谱
图，与原始信号对比，可以很明显的看出区别。
clear all;
clc;
fs=8000;
x1=wavread('pb8k.wav');
t=(0:length(x1)-1)/8000;
f=fs*(0:1023)/2048;
Au=0.05;
d=[Au*cos(2*pi*3800*t)]';%噪声为3.8kHz的余弦信号
x2=x1+d;
y1=fft(x1,2048);
y2=fft(x2,2048);
figure(1)
plot(t,x2)
grid on;axis tight;
title('加噪后的信号');
xlabel('time(s)');
ylabel('幅度');
figure(2)
subplot(2,1,1);
plot(f,abs(y1(1:1024)));grid on;axistight;
title('原始语音信号频谱');
xlabel('Hz');ylabel('幅度');
subplot(2,1,2);
plot(f,abs(y2(1:1024)));grid on;axistight;
title('加噪语音信号频谱');
xlabel('Hz');ylabel('幅度');
**********************************************************************************************************************
程序3:双线性变换法设计Butterworth滤波器
clear all;
clc;
fs=8000;
x1=wavread('pb8k.wav');
t=(0:length(x1)-1)/8000;
f=fs*(0:1023)/2048;
A1=0.05;A2=0.10;
d=[A1*cos(2*pi*3800*t)+A2*sin(2*pi*3600*t)]';
x2=x1+d;
wp=0.8*pi;
ws=0.85*pi;
Rp=1;
Rs=15;
Fs=8000;
Ts=1/Fs;
wp1=2/Ts*tan(wp/2);%将模拟指标转换成数字指标
ws1=2/Ts*tan(ws/2);
[N,Wn]=buttord(wp1,ws1,Rp,Rs,'s');%选择滤波器的最小阶数
[Z,P,K]=buttap(N);%创建butterworth模拟滤波器
[Bap,Aap]=zp2tf(Z,P,K);
[b,a]=lp2lp(Bap,Aap,Wn);
[bz,az]=bilinear(b,a,Fs);%用双线性变换法实现模拟滤波器到数字滤波器的转换
[H,W]=freqz(bz,az);%绘制频率响应曲线
figure(1)
plot(W*Fs/(2*pi),abs(H))
grid on;axis tight;
xlabel('频率（Hz）')
ylabel('频率响应')
title('Butterworth')
f1=filter(bz,az,x2);
figure(2)
subplot(2,1,1)
plot(t,x2)%画出滤波前的时域图
grid on;axis tight;
title('滤波前的时域波形');
subplot(2,1,2)
plot(t,f1);%画出滤波后的时域图
grid on;axis tight;
title('滤波后的时域波形');
y3=fft(f1,2048);
figure(3)
y2=fft(x2,2048);
subplot(2,1,1);
plot(f,abs(y2(1:1024)));%画出滤波前的频谱图
grid on;axis tight;
title('滤波前的频谱')
xlabel('Hz');
ylabel('幅度');
subplot(2,1,2)
plot(f,abs(y3(1:1024)));%画出滤波后的频谱图
grid on;axis tight;
title('滤波后的频谱')
xlabel('Hz');
ylabel('幅度');
%sound(x2,8000);
% sound(f1,8000);
********************************************************************************************************************
程序4:窗函数法设计滤波器：
clear all;
clc;
fs=8000;
x1=wavread('pb8k.wav');
t=(0:length(x1)-1)/8000;
f=fs*(0:2047)/4096;
A1=0.05;A2=0.10;
d=[A1*cos(2*pi*3600*t)+A2*sin(2*pi*3800*t)]';
x2=x1+d;
wp=0.8*pi;
ws=0.85*pi;
wdelta=ws-wp;
N=ceil(6.6*pi/wdelta);%取整
wn=(0.8+0.85)*pi/2;
[bz,az]=fir1(N,wn/pi,hamming(N+1));%选择窗函数，并归一化截止频率
figure(1)
freqz(bz,az);
grid on;axis tight;
f2=filter(bz,az,x2);
figure(2)
subplot(2,1,1)
plot(t,x2);
grid on;axis tight;
title('滤波前的时域波形');
subplot(2,1,2)
plot(t,f2);
grid on;axis tight;
title('滤波后的时域波形');
y3=fft(f2,4096);
f=fs*(0:2047)/4096;
figure(3)
y2=fft(x2,4096);
subplot(2,1,1);
plot(f,abs(y2(1:2048)));
grid on;axis tight;
title('滤波前的频谱')
xlabel('Hz');
ylabel('幅度');
subplot(2,1,2)
plot(f,abs(y3(1:2048)));
grid on;axis tight;
title('滤波后的频谱')
xlabel('Hz');
ylabel('幅度');
% sound(f2,8000);
***********************************************************************************************************************
程序5
clear all;close all;
t=0:0.01:3;
s=sin(2*pi*5*t)+sin(2*pi*15*t);
w=abs(fft(s));
w2=abs(fft(s,602));
figure;
subplot(211);plot(w);grid on; axistight;
subplot(212);plot(w2);grid on; axistight;
*****************************************************************************************************************************
程序6
clear all;close all;
X=[ones(1,51)zeros(1,60) ones(1,50)]; %Frequency response of an ideal low-pass filter
x=ifft(X); % IDFT of X withthe same number of points
Y=fft(x,483); % IDFT of Xafter it is padded with zeros
stem([1:161]/161-0.5,abs(fftshift(X)));% figure 1
axis([-0.5,0.5,0,1.5]);
figure;
stem(real(x)); % figure2
figure;
stem([1:483]/483-0.5,abs(fftshift(Y)));% figure 3
axis([-0.5,0.5,0,1.5]);
figure;
stem([1:161]/161-0.5,abs(fftshift(Y(1:3:481))));% figure 4
axis([-0.5,0.5,0,1.5]);
**********************************************************************************************************************
程序7
clear all;
clc
wp=0.8*pi;
ws=0.85*pi;
wdelta=ws-wp;
N=ceil(6.6*pi/wdelta);
wn=(0.8+0.85)*pi/2;
[bz,az]=fir1(N,wn/pi,hamming(N+1));
w=abs(fft(bz));w2=abs(fft(bz,399));
figure;
subplot(211);plot(w);grid on; axistight;
subplot(212);plot(w2);grid on; axistight;
-------------------------------------------------
h(1:201)=1;h(22:181)=0;
x1=ones(1,17);x2=ones(1,11);
w1=fft(x1,32);w2=fft(x2,32);
w3=w1.*w2; x3=ifft(w3);plot(x3);
