北京交通大学数字信号处理04DSP研究性学习报告多速率信号处理

1、数字信号处理课程研究性学习报告姓名 学号 同组成员 指导教师 时间 多速率信号处理专题研讨【目的】(1) 掌握序列抽取运算与内插运算的频谱变化规律。(2) 掌握确定抽取滤波器与内插滤波器的频率指标。(3) 掌握有理数倍抽样率转换的原理及方法。(4) 培养学生自主学习能力，以及发现问题、分析问题和解决问题的能力。【研讨题目】 基本题 1抽取、内插信号特征的时域/频域分析对于给定的单频模拟信号y(t)=sin(1000pt)，确定一个合适的采样频率fsam，获得离散信号yk，试进行以下问题的分析：(1) 对离散信号yk进行M=2倍抽取，对比分析yk和yMk在时域/频域的关系；(2) 对离散信号yk

2、进行L=2倍内插，对比分析yk和yk/L在时域/频域的关系。【温磬提示】在多速率信号分析中，离散序列的抽取和内插是多速率系统的基本运算，抽取运算将降低信号的抽样频率，内插运算将提高信号的抽样频率。两种运算的变换域描述中，抽取运算可能出现频谱线性混叠，而内插运算将出现镜像频谱。【设计步骤】1、 已知y(t)=sin(1000pt)频率为500Hz，周期为0.002s，可取时间范围T为0到0.004秒，两个周期，根据抽样定理取，每个周期抽取16个点。2、 用函数xD=x(1:M:end)对离散信号进行M=2倍的抽取，用fft计算频谱。3、 用函数xL=zeros(1,L*length(x);xL(

3、1:L:end)=x;对离散信号进行L=2的内插，用fft计算频谱。【仿真结果】对离散信号yk抽取和内插的时域/频域对比分析结果。抽取：内插：【结果分析】1、抽取运算在频域描述：对xk进行M倍抽取后得到的频谱为，即将xk的频谱扩展M倍，得到，再以为周期进行周期化并乘以因子。2、内插运算在频域描述：对xk进行L倍内插后得到频谱为，即将xk的频谱压缩L倍，除了与原序列频谱相差一个尺度因子外，两个序列频谱形状保持不变。由于的周期为，所以周期为，内插序列在的频谱将在区间重复L次。【发现问题】 (专题研讨或相关知识点学习中发现的问题)：应选择合适抽样频率来抽取足够多的点。【仿真程序】抽取：T=0.004

4、;fs=8000;Ts=1/fs;N=T/Ts;k=0:N-1;x=sin(1000*pi*(k.*Ts);subplot(2,2,1)stem(k.*Ts,x);xD=x(1:2:end);subplot(2,2,2);stem(xD);X=fft(x,32);X=fftshift(X);subplot(2,2,3);plot(abs(X);Y=fft(xD,16);Y=fftshift(Y);subplot(2,2,4);plot(abs(Y);内插：T=0.004;fs=8000;Ts=1/fs;N=T/Ts;k=0:N-1;x=sin(1000*pi*(k.*Ts);subplot(2

5、,2,1)stem(k.*Ts,x);xL=zeros(1,2*length(x);xL(1:2:end)=x;subplot(2,2,2);stem(xL);X=fft(x,32);X=fftshift(X);subplot(2,2,3);plot(abs(X);Y=fft(xL,64);Y=fftshift(Y);subplot(2,2,4);plot(abs(Y);【研讨题目】 中等题 2音乐信号多速率处理（1）分别给出抽样率为32KHz、16KHz和8KHz的三段音乐信号，利用Matlab仿真分别将信号通过抽样率为16KHz的D/A进行播放，试听播放结果分析其中存在的现象和问题；32K

6、Hz音乐信号文件：kdqg32k.wav16KHz音乐信号文件：kdqg16k.wav8 KHz音乐信号文件：kdqg8k.wav（2）设计多速率信号处理系统，使得抽样率为32KHz和8KHz的三段音乐信号通过抽样率为16KHz的D/A能够正常播放。 【设计步骤】（1） 用wavread函数读取文件，再使用sound函数对音乐信号以16Hz抽样频率播放。（2） 用wavread函数读取文件，对音乐信号进行内插或者抽取，再用sound函数播放。【仿真结果】（1）32kHz以16kHz播放声音低沉，16kHz能够正常播放，8kHz声音变得尖细。（2）音乐信号均正常播放。【结果分析】32kHz音乐信

7、号2倍抽取后，时域抽样间隔变为2倍，抽样频率变为原来二分之一，可以以16kHz正常播放；8kHz音乐信号2倍内插，时域抽样间隔变为二分之一，抽样频率变为原来2倍，可以以16kHz正常播放。【自主学习内容】wavread函数与sound函数，内插函数及抽取函数。【阅读文献】【发现问题】 (专题研讨或相关知识点学习中发现的问题)：【问题探究】【仿真程序】（1）x=wavread(C:UsersAdministratorDesktop数字信号研讨kdqg32k.wav)sound(x,16000);x=wavread(C:UsersAdministratorDesktop数字信号研讨kdqg16k.

8、wav)sound(x,16000);x=wavread(C:UsersAdministratorDesktop数字信号研讨kdqg8k.wav)sound(x,16000);（2）x=wavread(C:UsersAdministratorDesktop数字信号研讨kdqg32k.wav);y=decimate(x,2);sound(x,16000);x=wavread(C:UsersAdministratorDesktop数字信号研讨kdqg8k.wav);y=interp(x,2);sound(x,16000);【研讨题目】 实践题 3对于人体脉搏信号可以综合地反映人体的心脏器官和心血管

9、系统的生理变化。对人体脉搏信号进行非线性分析的研究，为评估人体心脏器官与血液循环系统的生理状态以及研究人的情绪对生理变化的影响提供了依据，并且对进一步深入认识人体生命系统的复杂性现象及规律具有重要意义。本题目将给出实际采集的人体脉搏信号，采样率为fsam=200Hz，疲劳状态下脉搏信号：疲劳.txt清晰状态下脉搏信号：清醒.txt(1) 利用Matlab读取脉搏信号，并进行脉搏信号的频谱分析，得出人体脉搏频谱特征。(2) 已知人体脉搏信号通常分布在1-20Hz内，设计合适的滤波器得到纯净的脉搏信号。(3) 由于脉搏信号通常分布在1-20Hz内，而人耳听觉范围在300-3400Hz范围，为使人耳

10、能够听到清醒的脉搏信号，请设计合适的多速率信号处理系统实现。【频谱特征分析】频谱主要部分是1-20Hz的低频分量。【设计步骤】1、 先画出脉搏信号的时域波形和频谱。2、 通过级联的巴特沃斯滤波器滤波3、 通过多速率信号处理系统实现对脉搏信号的播放【仿真结果】（1）（2）（3）【结果分析】1、 采用的脉搏信号中主要是低于20Hz的低频分量，在滤波是采用巴特沃斯高通滤波器和巴特沃斯低通滤波器的级联，比带通滤波器实现的代价要低。2、 通过内插滤波器可以实现对低频率信号的播放，内插产生镜像分量，取出可听的频率分量。【仿真程序】（1）fs=200;T=1/fs;x1=load(C:UsersAdmini

11、stratorDesktop数字信号研讨清醒.txt);x2=load(C:UsersAdministratorDesktop数字信号研讨疲劳.txt);M1=length(x1);M2=length(x2);t1=(0:M1-1)*T;t2=(0:M2-1)*T;N1=M1*2-1;N2=M2*2-1;X1=fftshift(fft(x1,N1);X2=fftshift(fft(x2,N2);W1=(-M1+1:M1-1)*fs/N1;W2=(-M2+1:M2-1)*fs/N2;subplot(2,2,1);plot(t1,x1);title(清醒状态脉搏信号时域波形);xlabel(t(s

12、);subplot(2,2,2);plot(t2,x2);title(疲劳状态脉搏信号时域波形);xlabel(t(s);subplot(2,2,3);plot(W1,abs(X1);title(清醒状态脉搏信号频谱);xlabel(频率);ylabel(幅度);axis(0,100,0,3.5*105);subplot(2,2,4);plot(W2,abs(X2);title(疲劳状态脉搏信号频谱);xlabel(频率);ylabel(幅度);axis(0,100,0,3.5*105);（2）fs=200;T=1/fs;x=load(C:UsersAdministratorDesktop数字

13、信号研讨清醒.txt);M=length(x);t=(0:M-1)*T;N=M*2-1;W=(-M+1:M-1)*fs/N;%Butterworth高通滤波器Wp1=0.9/(fs/2);Ws1=0.1/(fs/2);n,Wn=buttord(Wp1,Ws1,0.5,40);b,a=butter(n,Wn,high);y=filtfilt(b,a,x);%Butterworth低通滤波器Wp=30/(fs/2);Ws=50/(fs/2);Rp=3;Rs=60;n,Wc=buttord(Wp,Ws,Rp,Rs);b,a=butter(n,Wc,low);q=filter(b,a,y);Y=fft

14、shift(fft(q,N);figure(1);plot(t,q);title(滤波后清醒状态脉搏信号时域波形);xlabel(t(s);figure(2);plot(W,abs(Y);title(滤波后清醒状态脉搏信号频谱);xlabel(频率);ylabel(幅度);axis(0,100,0,15000);fs=200;T=1/fs;x=load(C:UsersAdministratorDesktop数字信号研讨疲劳.txt);M=length(x);t=(0:M-1)*T;N=M*2-1;W=(-M+1:M-1)*fs/N;%Butterworth高通滤波器Wp1=0.9/(fs/2)

15、;Ws1=0.1/(fs/2);n,Wn=buttord(Wp1,Ws1,0.5,40);b,a=butter(n,Wn,high);y=filtfilt(b,a,x);%Butterworth低通滤波器Wp=30/(fs/2);Ws=50/(fs/2);Rp=3;Rs=60;n,Wc=buttord(Wp,Ws,Rp,Rs);b,a=butter(n,Wc,low);q=filter(b,a,y);Y=fftshift(fft(q,N);figure(1);plot(t,q);title(滤波后清醒状态脉搏信号时域波形);xlabel(t(s);figure(2);plot(W,abs(Y)

16、;title(滤波后清醒状态脉搏信号频谱);xlabel(频率);ylabel(幅度);axis(0,100,0,15000);（3）%以下通过内插实现对脉搏信号的播放%先对人清醒状态的脉搏信号进行带通滤波，滤掉噪声信号fs=200;T=1/fs;x=load(C:UsersAdministratorDesktop数字信号研讨清醒.txt);%Butterworth高通滤波器Wp1=0.9/(fs/2);Ws1=0.1/(fs/2);n,Wn=buttord(Wp1,Ws1,0.5,40);b,a=butter(n,Wn,high);y=filtfilt(b,a,x);%Butterworth

17、低通滤波器Wp=30/(fs/2);Ws=50/(fs/2);Rp=3;Rs=60;n,Wc=buttord(Wp,Ws,Rp,Rs);b,a=butter(n,Wc,low);q=filter(b,a,y);W=linspace(0,pi,1024);Q=freqz(q,1,W);figure(1);plot(W/pi,abs(Q);title(原始脉搏信号去噪后的频谱);xlabel(归一化频率);ylabel(幅度);%对滤波后的噪音进行L倍内插L=25;y1=upsample(q,L);Y1=freqz(y1,1,W);figure(2);plot(W/pi,abs(Y1);title(25倍内插后脉搏信号的频谱);xlabel(归一化频率);ylabel(幅度);%构件带