数字信号处理课程设计程序语音信号处理
%请将此文件设置为.m文件
%数字信号处理课程设计程序如下：
%原始语音信号的波形和频谱：
[x1,fs,nbits]=wavread ('215');%fs为采样频率
N=length (x1) ; %求出语音信号的长度
Y=fft(x1,N); %傅里叶变换
sound(x1,fs,nbits); %回放语音信号
subplot(2,1,2);%两行两列的第二个图形
plot(abs(Y));subplot(2,1,1);%abs表示取绝对值，两行两列的第一个图形
plot(x1);title('原始信号波形');%plot显示括号里面的图形
title('原始信号频谱')
%FFT（快速傅里叶）变换：
y1=fft(x1,1024);%对采样后信号做1024点FFT变换
f=fs*(0:511)/1024;%生成频率序列
figure(1)%表示合在一个图里面
subplot(2,1,1);%两行一列的第一个图形
plot(abs(y1(1:512)));%原始语音信号的频谱
title('原始语音信号频谱')
xlabel('频率/Hz');
ylabel('幅值')
subplot(2,1,2);%两行一列的第二个图
plot(abs(y1(1:1024))) %采样后信号的FFT频谱图
title('原始语音信号FFT频谱')
xlabel('点数N');
ylabel('幅值');
%双线性变换法设计IIR数字低通滤波器：
fs=44100;
[x2]=wavread('215');
Ts=1/fs;
wp=2*pi*1000/fs; %通带截止频率
ws=2*pi*1200/fs; %阻带截止频率
Rp=1; %通带衰减
Rl=10; %阻带衰减
wp1=2/Ts*tan(wp/2);%课本P187页，将数字指标转换成模拟指标
ws1=2/Ts*tan(ws/2);
[N,Wn]=buttord(wp1,ws1,Rp,Rl,'s'); %选择滤波器的最小阶数
[Z,P,K]=buttap(N);
[Bap,Aap]=zp2tf(Z,P,K);%创建butterworth模拟低通滤波器
[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
xlabel('频率／Hz')
ylabel('频率响应幅度')
title('IIR低通滤波器')
f1=filter(bz,az,x2);%进行滤波处理
figure(2)
subplot(2,1,1)
plot(x2) %画出滤波前的时域图
title('IIR低通滤波器滤波前的时域波形');
subplot(2,1,2)
plot(f1); %画出滤波后的时域图
title('IIR低通滤波器滤波后的时域波形');
sound(f1,44100); %播放滤波后的信号
F0=fft(f1,1024);
f=fs*(0:511)/1024;
figure(3)
y2=fft(x2,1024);
subplot(2,1,1);
plot(f,abs(y2(1:512))); %画出滤波前的频谱图
title('IIR低通滤波器滤波前的频谱')
xlabel('频率/Hz');
ylabel('幅值');
subplot(2,1,2)
F1=plot(f,abs(F0(1:512))); %画出滤波后的频谱图
title('IIR低通滤波器滤波后的频谱')
xlabel('频率/Hz');
ylabel('幅值');
%窗函数法设计FIR数字低通滤波器：
fs=44100;
wp=2*pi*1000/fs;
ws=2*pi*1200/fs;
Rp=1;
Rs=100;
wdelta=ws-wp;
N=ceil(8*pi/wdelta); %取整
wn=(wp+ws)/2;
[b,a]=fir1(N,wn/pi,hamming(N+1)); %选择窗函数，并归一化
截止频率
figure(1)
freqz(b,a,512);%绘制频率响应函数
title('FIR低通滤波器');
f2=filter(b,a,x1);%进行滤波处理
figure(2)
subplot(2,1,1)
plot(x1)
title('FIR低通滤波器滤波前的时域波形');
subplot(2,1,2)
plot(f2);
title('FIR低通滤波器滤波后的时域波形');
sound(f2,44100); %播放滤波后的语音信号
F0=fft(f2,1024);
f=fs*(0:511)/1024;
figure(3)
y2=fft(x1,1024);
subplot(2,1,1);
plot(f,abs(y2(1:512)));
title('FIR低通滤波器滤波前的频谱')
xlabel('频率/Hz');
ylabel('幅值');
subplot(2,1,2)
%双线性设计高通滤波器cheby1（切比雪夫滤波器）：
plot(x1);
title('采样后语音信号时域波形');%？？？
y1=fft(x1);
figure;
plot(abs(y1));
title('采样后语音信号频谱图');
fp=5000;fc=4800;As=100;Ap=1;Fs=44100;
wc=2*pi*fc/Fs;wb=2*pi*fp/Fs;%保证在0-1之间
[N,wp]=cheb1ord(wb,wc,Ap,As);
[b,a]=cheby1(N,Ap,wp,'high');
figure;
freqz(b,a);%绘制频率响应曲线
x=filter(b,a,x1);%进行滤波处理
X=fft(x,8192);%？？？
figure;
subplot (2,2,1);plot (abs(y1));%y1为x1的fft变换
title('滤波前信号频谱');
subplot (2,2,2);plot (abs(X));
title('滤波后信号频谱');
subplot(2,2,3);plot(x1);
title('滤波前信号波形');
subplot(2,2,4);plot(x);
title ('滤波后信号波形');
sound(x1,44100,16);
%窗函数法设计FIR数字滤波器（高通）：
fs=44100;
wp=2*pi*4800/fs;
ws=2*pi*5000/fs;
Rp=1;
Rs=100;
wdelta=ws-wp;
N=ceil(8*pi/wdelta); %取整
wn=(wp+ws)/2;
[b,a]=fir1(N,wn/pi,hamming(N+1)); %选择窗函数，并归一化截止频率
figure(1)
freqz(b,a,512);%绘制频率响应图形
title('FIR高通滤波器');
f2=filter(b,a,x1);%进行滤波处理
figure(2)
subplot(2,1,1)
plot(x1)
title('FIR高通滤波器滤波前的时域波形');
subplot(2,1,2)
plot(f2);
title('FIR高通滤波器滤波后的时域波形');
sound(f2,44100); %播放滤波后的语音信号
F0=fft(f2,1024);%对f2进行1024点fft变换
f=fs*(0:511)/1024;
figure(3)
y2=fft(x1,1024);%对x1进行1024点fft变换
subplot(2,1,1);
plot(f,abs(y2(1:512)));
title('FIR高通滤波器滤波前的频谱')
xlabel('频率/Hz');
ylabel('幅值');
subplot(2,1,2)
F2=plot(f,abs(F0(1:512)));
title('FIR高通滤波器滤波后的频谱')
xlabel('频率/Hz');
ylabel('幅值');
%双线性设计带通滤波器ellip：
plot(x1);
title('采样后语音信号时域波形');%？？？
y1=fft(x1);%对x1
figure;
plot(abs(y1));
title('采样后语音信号频谱图');
fp1=1200;fp2=3000;fc1=1000;fc2=3200;As=100;Ap=1;Fs=8000;
wc=[2*fc1/Fs,2*fc2/Fs];wp=[2*fp1/Fs,2*fp2/Fs];
[N,Wn]=ellipord (wp,wc,Ap,As);
[b,a]=ellip(N,Ap,As,Wn,'stop');
[h,f]=freqz(b,a,2^6,Fs);
figure;
plot(f,abs(h));
title('滤波器幅频响应图');
p=angle(h);%相频响应
ph=p*180/pi;%？？？
figure;
plot(f,ph);
title('滤波器相频响应图');
y2=filter(b,a,x1);
figure;
plot(y2);
t
itle('滤波后语音信号时域波形');
y3=fft(y2);
figure ;
plot(abs(y3));
title('滤波后语音信号频谱图');
sound(y2,44100);
%窗函数设计带通滤波器：
y1=fft(x1);
fp1=1200;fp2=3000;fc1=1000;fc2=3200;As=100;Ap=1;Fs=8000;
wp1=2*pi*fp1/Fs;wc1=2*pi*fc1/Fs;wp2=2*pi*fp2/Fs;wc2=2*pi*fc2/Fs;
wdel=wp1-wc1;
beta=0.112*(As-8.7);%？？？
N=ceil ((As-8)/2.285/wdel);%？？？
wn=kaiser(N+1,beta);%？？？
ws=[(wp1+wc1)/2/pi,(wp2+wc2)/2/pi];
b=fir1(N,ws,wn);
figure;
freqz(b,1);%绘制频率响应
x=fftfilt (b,x1);%？？？
X=fft(x,8192);%？？？
figure;
subplot (2,2,1);plot (abs(y1));
title('滤波前信号频谱');
subplot (2,2,2);plot (abs(X));
title('滤波后信号频谱');
subplot(2,2,3);plot(x1);
title('滤波前信号波形');
subplot(2,2,4);plot(x);
title ('滤波后信号波形');
sound(x1,44100,16);