04北交大数字信号处理多速率

《数字信号处理》课程研究性学习报告姓名张然学号13211074同组成员蔡逸飞13211078朱斌指导教师陈后金时间2015/6多速率信号处理专题研讨【目的】(1)掌握序列抽取运算与内插运算的频谱变化规律。(2)掌握确定抽取滤波器与内插滤波器的频率指标。(3)掌握有理数倍抽样率转换的原理及方法。(4)培养学生自主学习能力，以及发现问题、分析问题和解决问题的能力。【研讨题目】基本题1．抽取、内插信号特征的时域/频域分析对于给定的单频模拟信号y(t)=sin(1000t)，确定一个合适的采样频率fsam，获得离散信号y[k]，试进行以下问题的分析：(1)对离散信号y[k]进行M=2倍抽取，对比分析y[k]和y[Mk]在时域/频域的关系；(2)对离散信号y[k]进行L=2倍内插，对比分析y[k]和y[k/L]在时域/频域的关系。【温磬提示】在多速率信号分析中，离散序列的抽取和内插是多速率系统的基本运算，抽取运算将降低信号的抽样频率，内插运算将提高信号的抽样频率。两种运算的变换域描述中，抽取运算可能出现频谱混叠线性，而内插运算将出现镜像频谱。【设计步骤】1.对所给定信号进行抽样，由于所给定信号为y(t)=sin(1000t)，其频率f=500Hz,且实验过程中我们发现，当T取较大值时，Matlab绘出的图像会有较大误差，经查询发现这是由于Matlab的数值计算特点所决定的，于是我们减小了T的值，同时加大了采样频率fsam,，使其能产生理想的图像。2.对离散信号进行M=2倍的抽取时，我们通过xD=x(1:M:end);进行MATLAB的仿真运算。3.对离散信号进行L=2倍内插时，我们通过xL=zeros(1,L*length(x));xL(1:L:end)=x;进行MATLAB的仿真运算。4.对于各离散信号的频谱计算，我们都采用了FFT进行计算。【仿真结果】对离散信号y[k]抽取和内插的时域/频域对比分析结果。【结果分析】抽取、内插运算在频域的描述和性质。通过图像我们可以发现，在时域内，离散序列的抽取表示每隔M-1点抽取原序列中的一点；离散信号的内插是在原序列每两个样本点之间插入L-1个零值样本点。在频域内，M倍抽取后序列的频谱为将原序列频谱扩展M倍，得到，再将其周期化后的函数相加并乘以因子1/M，即可得到抽取后序列的频谱。由于原序列频谱在区间[-π，π]范围内满足，所以抽取后信号的频谱不会发生混叠。信号L倍内插的频谱为将原序列频谱压缩L倍得到，除了与原序列频谱相差一个尺度因子外，两个频谱的形状保持不变。由于的周期为2π，所以的周期为2π/L。内插序列在区间[-π/L，π/L]内的频谱将在区间[-π，π]内重复L次，此重复的部分为镜像频谱【自主学习内容】【阅读文献】[1]陈后金.数字信号处理[M].北京：高等教育出版社.2008.11【发现问题】(专题研讨或相关知识点学习中发现的问题)：当连续信号频率过大时，T值较大时MATLAB绘图会产生较大误差【问题探究】【仿真程序】抽取：T=0.004;fs=8000;Ts=1/fs;N=T/Ts;x=zeros(1,N);t=0:N-1;x=sin(1000*pi*(t*Ts));subplot(2,2,1);stem(t*Ts,x);xD=x(1:2:end);subplot(2,2,3);stem(xD);X=fft(x,32);X=fftshift(X);subplot(2,2,2);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;x=zeros(1,N);t=0:N-1;x=sin(1000*pi*(t*Ts));subplot(2,2,1);stem(t*Ts,x);xL=zeros(1,2*length(x));xL(1:2:end)=x;subplot(2,2,3);stem(xL);X=fft(x,32);X=fftshift(X);subplot(2,2,2);plot(abs(X));Y=fft(xD,64);Y=fftshift(Y);subplot(2,2,4);plot(abs(Y));【研讨题目】中等题2．音乐信号多速率处理（1）分别给出抽样率为32KHz、16KHz和8KHz的三段音乐信号，利用Matlab仿真分别将信号通过抽样率为16KHz的D/A进行播放，试听播放结果分析其中存在的现象和问题；32KHz音乐信号文件：kdqg32k.wav16KHz音乐信号文件：kdqg16k.wav8KHz音乐信号文件：kdqg8k.wav（2）设计多速率信号处理系统，使得抽样率为32KHz和8KHz的三段音乐信号通过抽样率为16KHz的D/A能够正常播放。D/AA/Dfsam=32kHzx(t)x[k]y(t)fsam=16kHzH(z)y[k]D/AA/Dfsam=8kHzx(t)x[k]y(t)fsam=16kHzH(z)y[k]【设计步骤】(1)用wavread函数读取文件，再使用sound函数用16k对信号进行恢复播放(2)用wavread函数读取文件，在对信号用decimate或interp函数变换频率，再用sound函数播放。【仿真结果】（1）第一个8kHz的音频用16k播放的时候声音变得尖细。16k用16k播放时声音正常，32k用16k播放时声音变得粗狂。（2）进行内插或抽取后的信号播放均正常。【结果分析】8kHz的文件，在内插之后，时域展宽。对应的频域减小，fs增大，再用16k恢复使信号正常播放。32kHz刚好相反。因此对速率处理系统可以很方便的实现速率转换。【自主学习内容】Wavread函数与sound函数的使用，抽取内插的改变播放速率的原理【阅读文献】【发现问题】(专题研讨或相关知识点学习中发现的问题)：【问题探究】【仿真程序】（1）y=wavread('C:\Users\Administrator\Desktop\kdqg8k.wav');sound(y,16000);y=wavread('C:\Users\Administrator\Desktop\kdqg16k.wav');sound(y,16000);y=wavread('C:\Users\Administrator\Desktop\kdqg32k.wav');sound(y,16000);（2）y=wavread('C:\Users\Administrator\Desktop\kdqg32k.wav');x=decimate(y,2);sound(x,16000)y=wavread('C:\Users\Administrator\Desktop\kdqg8k.wav');x=interp(y,2);sound(x,16000)【研讨题目】实践题3．对于人体脉搏信号可以综合地反映人体的心脏器官和心血管系统的生理变化。对人体脉搏信号进行非线性分析的研究，为评估人体心脏器官与血液循环系统的生理状态以及研究人的情绪对生理变化的影响提供了依据，并且对进一步深入认识人体生命系统的复杂性现象及规律具有重要意义。本题目将给出实际采集的人体脉搏信号，采样率为fsam=200Hz，疲劳状态下脉搏信号：疲劳.txt清晰状态下脉搏信号：清醒.txt(1)利用Matlab读取脉搏信号，并进行脉搏信号的频谱分析，得出人体脉搏频谱特征。(2)已知人体脉搏信号通常分布在1-20Hz内，设计合适的滤波器得到纯净的脉搏信号。(3)由于脉搏信号通常分布在1-20Hz内，而人耳听觉范围在300-3400Hz范围，为使人耳能够听到清醒的脉搏信号，请设计合适的多速率信号处理系统实现。【频谱特征分析】频谱的主要部分是1-20Hz的低频分量【设计步骤】1.先求出脉搏信号的时域波形和频谱2.通过级联的巴特沃斯滤波器滤波3.通过多速率信号处理系统实现对脉搏信号的播放【仿真结果】【结果分析】1.采集的脉搏信号中主要是低于20Hz的低频分量，在滤波是采用巴特沃斯高通滤波器和波特沃斯低通滤波器的级联，比带通滤波器实现的代价要低2.滤波的效果还算理想，但由于有过渡带的存在，高频部分滤的不干净3.通过内插滤波器可以实现对低频率信号的播放，内插产生的镜像分量，取出可听的频率分量【仿真程序】%以下对人清醒状态/疲劳状态的脉搏信号进行时/频分析fs=200;T=1/fs;x1=load('lucid.txt');%导入清醒状态的脉搏信号x2=load('fatigue.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)');subplot(2,2,2);plot(t2,x2);title('疲劳状态脉搏信号时域波形');xlabel('t(s)');%输出脉搏信号的频谱subplot(2,2,3);plot(W1,abs(X1));title('清醒状态脉搏信号频谱');xlabel('Frequency(Hz)');ylabel('Amplitude');axis([0,100,0,3.5*10^5]);subplot(2,2,4);plot(W2,abs(X2))title('疲劳状态脉搏信号频谱');xlabel('Frequency(Hz)');ylabel('Amplitude');axis([0,100,0,3.5*10^5]);%以下对人清醒状态的脉搏信号进行滤波fs=200;T=1/fs;x=load('lucid.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=fftshift(fft(q,N));figure(1);plot(t,q);title('滤波后清醒状态脉搏信号时域波形');xlabel('t(s)');figure(2);plot(W,abs(Y));title('滤波后清醒状态脉搏信号频谱');xlabel('Frequency(Hz)');ylabel('Amplitude');axis([0,100,0,15000]);%以下对人疲劳状态的脉搏信号进行滤波fs=200;T=1/fs;x=load('fatigue.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,W