实验5：MATLAB在数字信号处理中的应用(4学时）
掌握MATLAB对离散时间信号与系统的时域和频域描述及分析
【实验目的】
1、 掌握离散时间信号与系统的时域和频域
分析方法
【实验内容】
1、 练习书上实例，验证结果
2、录制一段各人自己的语音信号,存成.wav文件,录制时间不要太长，5-10秒即可。利用MATLAB提供的wavread函数进行读取，而后对录制的信号x(t)进行采样；画出采样后语音信号的时域波形和频谱图。（提示，仿例题7.10,选择不同的采样频率fs,观察频谱）
3、仿照教材7.1节对该语音信号进行各种时域处理（如截取，合成，添加随机噪声，周期延拓，移位），观察波形图，若可能可将其用MATLAB存成wav文件，听一下。（提示：函数randn可构造随机噪声序列；
人的语音信号的范围20Hz-20KHz,选择合适的采样频率，在数字通信中，根据采样定理, 最小采样频率为语音信号最高频率的2倍。电话通令的频带是300-3400Hz，所以最低采样频率为6800Hz。根据CCITT的长途通信协议的标准，目前各国一般采用的声音采样频率为8kHz；所以，为了将来系统与其他通信设备交流的方便，最终采样频率一般确定为8000Hz)。
4、用FFT函数对上述语音信号，进行频域分析，观察该语音信号的频谱范围.（仿例题7.19）.
5、构造一个理想数字滤波器，要求保留1000HZ-2000HZ以内的语音信号，观察信号波形图。 【实验报告】 Q701
x1=[0,1,2,3,4,3,2,1,0];ns1=-2; x2=[2,2,0,0,0,-2,-2]; ns2=2; nf1=ns1+length(x1)-1; nf2=ns2+length(x2)-1;
ny= min(ns1,ns2):max(nf1,nf2); a1 = zeros(1,length(ny)); xa2 = xa1;
xa1(find((ny>=ns1)&(ny<=nf1)==1))=x1;
xa2(find((ny>=ns2)&(ny<=nf2)==1))=x2;
ya = xa1 + xa2 yp = xa1.* xa2 subplot(4,1,1), stem(ny,xa1,'.') subplot(4,1,2), stem(ny,xa2,'.') line([ny(1),ny(end)],[0,0]) subplot(4,1,3), stem(ny,ya,'.') line([ny(1),ny(end)],[0,0]) subplot(4,1,4), stem(ny,yp,'.') line([ny(1),ny(end)],[0,0])
clear;close all
n0=input('输入序列起点:n0='); N=input('输入序列长度:N='); n1=input('输入位移:n1='); n=n0:n1+N+5; u=[(n-n1)>=0];
x1=[(n-n1)>=0]-[(n-n1-N)>=0]
x2=[(n>=n1)&(n<(N+n1))]; x3=exp(j*n*pi/8).*x2;
subplot(2,2,1);stem(n,x1,'.');
xlabel('n');ylabel('x1(n)');
axis([n0,max(n),0,1]);
subplot(2,2,3);stem(n,x2,'.');
xlabel('n');ylabel('x2(n)');
axis([n0,max(n),0,1]);
subplot(2,2,2);stem(n,real(x3),'.');
xlabel('n');ylabel('x3(n)的实部');
line([n0,max(n)],[0,0]);
axis([n0,max(n),-1,1]);
subplot(2,2,4);stem(n,imag(x3),'.');
xlabel('n');ylabel('x3(n)的虚部');
line([n0,max(n)],[0,0]);
axis([n0,max(n),-1,1]);
输入序列起点:n0=-6
输入序列长度:N=15
输入位移:n1=3
x1 = Columns 1 through 13 0 0 0 0 0 0 0 0 0 1 1 1 1
Columns 14 through 26
1 1 1 1 1 1 1 1 1 1 1 0 0
Columns 27 through 30
0 0 0 0
Q03
clear;close all N=24;M=8;
m=input('输入移位值：m='); if (m<1|m>=N-M+1)
fprintf('输入数据不在规定范围内！');
break end
n=0:N-1;
x1=(0.8).^n; x2=[(n>=0)&(n xn=x1.*x2; xm=zeros(1,N); for k=m+1:m+M xm(k)=xn(k-m); end xc=xn(mod(n,8)+1); xcm=xn(mod(n-m,8)+1); subplot(3,2,1); stem(n,xn,'.'); title('x(n)');axis([0,length(n),0,1]); ylabel('x(n)');xlabel('n'); subplot(3,2,2);stem(n,xm,'.'); title('x(n-m)');axis([0,length(n),0,1]); ylabel('x(n-m)');xlabel('n'); subplot(3,2,5); stem(n,xc,'.'); title('x(n)的周期延拓');axis([0,length(n),0,1]); ylabel('x(mod(n,8))');xlabel('n'); subplot(3,2,6);stem(n,xcm,'.'); title('x(n)的循环移位');axis([0,length(n),0,1]); ylabel('x(mod(n-m,8))');xlabel('n'); 输入移位值：m=2 试验一 语音信号处理 语音信号处理综合运用了数字信号处理的理论知识，对信号进行计算及频谱分析，设计滤波器，并对含噪信号进行滤波。 一，具体分为以下步骤： （1）语音信号的采集：利用Windows 下的录音机，录制一段话音。然后在Matlab软件平台下，利用函数wavread对语音信号进行采样，播放语音信号，并绘制原始语音信号； （2）对原始信号加入噪声：对原始语 音信号加入 s=0.05*sin(2*pi*f*Ts*n) 的噪声，采样后可知Fs ，选择f = 2500， 播放加入噪声信号的语音信号，并绘制噪声信号和含噪语音信号； （3）频谱分析：分别对原始语音信号， 噪声信号和含噪声的语音信号进行频谱分析，并绘出各频谱图； （4）设计滤波器：计算滤波器的性能 指标，设计滤波器，绘制滤波器的特性曲线； （5）滤波器滤波：用自己设计的滤波 器对采集的信号进行滤波，得出滤波后信号 的时域波形和频谱，并对滤波前后的信号进行对比，分析信号的变化，并回放语音信号，感觉滤波前后的声音有变化。 (6) 对原始信号进行整数倍抽取，比较 抽取前后的频谱图 （7）对原始信号进行整数倍内插，比较原始信号频谱，内插零值时的频谱和滤波后的频[x,fs,bits]=wavread('aaa.wav'); N=length(x); f=2500; x=x'; n=0:1/fs:(N-1)/fs; s=0.05*sin(2*pi*f*n); y=x+s; f1=0:fs/N:fs*(N-1)/N; x1=fft(x,N); s1=fft(s,N); Y=fft(y,N); figure(1); subplot(311);plot(n,x);title('原始信号波形图'); xlabel('t/s');ylabel('幅值'); subplot(312);plot(n,s);title('噪声信号波形图');ylim([-0.08,0.08]); xlabel('t/s');ylabel('幅值'); subplot(313);plot(n,y);title('信号+噪声信号波形图'); xlabel('t/s');ylabel('幅值'); figure(2); subplot(221);plot(f1,abs(x1));title('原始信号频谱图'); xlabel('频率/Hz'); ylabel('频率响应幅度');xlim([0,fs/12]); subplot(222);plot(f1,abs(s1));title('噪声信号频谱图'); xlabel('频率/Hz'); ylabel('频率响应幅度');xlim([0,fs/12]); subplot(223);plot(f1,abs(Y));title('信号+噪声频谱图'); xlabel('频率/Hz'); ylabel('频率响应幅度 ');xlim([0,fs/12]); 谱图。 滤波器设计 wp=1500*2*pi;ws=2200*2*pi;Rp=2;Rs=15; [N1,Wn]=buttord(wp,ws,Rp,Rs,'s'); %创建butterworth模拟滤波器 [Z,P,K]=buttap(N1); %把滤波器的零极点模型转化为传递函数模型 [Bap,Aap]=zp2tf(Z,P,K); [b,a]=lp2lp(Bap,Aap,Wn); [bz,az]=bilinear(b,a,fs); [H,W]=freqz(bz,az); figure(3); subplot(211), plot(W*fs/(2*pi),abs(H)); grid title('butterworth低通滤波器频响图');xlabel('频率/Hz'); ylabel('频率响应幅度'); y1=filter(bz,az,y); subplot(212),plot(n,y1);title('滤波后信号波形 图'); xlabel('频率/Hz'); ylabel('频率响应幅度'); 原始信号波形图0.2值幅0-0.200.10.20.30.40.5t/s噪声信号波形图0.05值幅0-0.0500.10.20.30.40.5t/s信号+噪声信号波形图0.2值幅0-0.200.10.20.30.40.5t/s 对原始信号进行整数倍抽取 M=2; for i=1:fix(N/M) y2(i)=x(M*i); end y2=fft(y2,i); fs2=fs/M; f2=0:fs2/i:fs2*(i-1)/i; figure(4);subplot(322);plot(f2,abs(y2)); title('原始信号');xlim([0 fs/(12*M/2^1)]); M1=1;%%按M1倍内插 y3=zeros(1,N*M1); for i=1:N y3(i*M1)=x(i); end y3=fft(y3,N*M1);%%内插后进行傅里叶变换 fs3=fs*M1; f3=0:fs3/(N*M1):fs3*(N*M1-1)/(N*M1); figure(5);subplot(221);plot(f3,abs(y3)); title('原始信号频谱图');xlim([0 150000]); 原始信号频谱图噪声信号频谱图30015000.6度0.7200度幅幅1000应应响100响率率500频频002000400000200040000.60.7频率/Hz频率/Hz信号+噪声频谱图1500度幅1000应响率500频0.600.7020004000频 率/Hzbutterworth低通滤波器频响图1.5度幅1应响率0.5频000.511.522.533.5频率/Hzx 104滤波后信号波形图0.2度0.1幅应响0率频-0.1-0.200.10.20.30.40.50.60.7频率/Hz按2倍抽取后的频谱按4倍抽取后的频谱2001001005000200040000020004000按8倍抽取后的频谱按16倍抽取后的频谱4020201000200040000020004000按32倍抽取后的频谱按64倍抽取后的频谱1045200200040000020004000 四，结果分析 对原始信号加入噪声后播放，声音中明显存在单频尖锐噪声。后对信号+噪声信号进行滤波后，播放的声音中仍留有少量的噪声信号，且整体信号不如加噪前明亮，原因是采用低通滤波器滤去高频噪声时，原信号的高频信号也被滤去了，所以声音没加噪前丰富，显得低沉。 对原始信号进行整数倍抽取，相当于降低采样频率，由图可知，当抽取倍数为16倍时出现了明显的频谱混叠现象。原因是由于抽样频率fs不满足fs>=2*fh，从而周期延拓时出现了频谱混叠现象。而对信号进行整数倍内插则没有使信号频谱变差。从而也验证了采样定理的正确性。 原始信号频谱图30020010003002001000按2倍内插零值后的频谱05101540510154x 10按4倍内插零值后的频谱3002001000x 10按8倍内插零值后的频谱3002001000051015x 104051015x 104