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

《数字信号处理》课程研究性学习报告

姓名

学号

同组成员

指导教师

时间

多速率信号处理专题研讨

【目的】

(1) 掌握序列抽取运算与内插运算的频谱变化规律。 (2) 掌握确定抽取滤波器与内插滤波器的频率指标。 (3) 掌握有理数倍抽样率转换的原理及方法。

(4) 培养学生自主学习能力，以及发现问题、分析问题和解决问题的能力。 【研讨题目】 基本题

1．抽取、内插信号特征的时域/频域分析

对于给定的单频模拟信号y (t )=sin(1000πt )，确定一个合适的采样频率f sam ，获得离散信号y [k ]，试进行以下问题的分析：

(1) 对离散信号y [k ]进行M=2倍抽取，对比分析y [k ]和y [M k ]在时域/频域的关系； (2) 对离散信号y [k ]进行L=2倍内插，对比分析y [k ]和y [k /L]在时域/频域的关系。

【温磬提示】

在多速率信号分析中，离散序列的抽取和内插是多速率系统的基本运算，抽取运算将降低信号的抽样频率，内插运算将提高信号的抽样频率。两种运算的变换域描述中，抽取运算可能出现频谱线性混叠，而内插运算将出现镜像频谱。 【设计步骤】

1、 已知y (t )=sin(1000πt )频率为500Hz ，周期为0.002s ，可取时间范围T 为0到0.004秒，两个周期，

根据抽样定理取Hz f sam 8000=，每个周期抽取16个点。

2、 用函数xD=x(1:M:end)对离散信号进行M=2倍的抽取，用fft 计算频谱。

3、 用函数xL=zeros(1,L*length(x));xL(1:L:end)=x;对离散信号进行L=2的内插，用fft 计算频谱。 【仿真结果】

对离散信号y [k ]抽取和内插的时域/频域对比分析结果。 抽取：

内插：

【结果分析】

1、抽取运算在频域描述：对x[k]进行M 倍抽取后得到][k x D 的频谱为∑-=-ΩΩ

=

10

2)(1

)(M l M

l j

j D e

X M

e X π，

即将x[k]的频谱)(Ω

j e

X 扩展M 倍，得到)(M

j e X Ω

，再以π2为周期进行周期化并乘以因子M 1。

2、内插运算在频域描述：对x[k]进行L 倍内插后得到][k x L 频谱为)()(ΩΩ=jL j L e X e X ，即将x[k]的频谱)(Ωj e X 压缩L 倍，除了与原序列频谱相差一个尺度因子外，两个序列频谱形状保持不变。由于)(Ωj e X 的周期为π2，所以)(ΩjL e X 周期为L

π

2，内插序列在],[L L ππ-的频谱将在区间

],[ππ-重复L 次。

【发现问题】 (专题研讨或相关知识点学习中发现的问题)： 应选择合适抽样频率来抽取足够多的点。 【仿真程序】

抽取： T=0.004; 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,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 进行播放，试听播放结果分析其中存在的现象和问题；

32KHz 音乐信号文件：kdqg32k.wav 16KHz 音乐信号文件：kdqg16k.wav 8 KHz 音乐信号文件：kdqg8k.wav （2）设计多速率信号处理系统，使得抽样率为32KHz 和8KHz 的三段音乐信号通过抽样率为16KHz 的D/A 能够正常播放。

sam x (t )sam

sam x (t )

sam

【设计步骤】

（1）用wavread 函数读取文件，再使用sound 函数对音乐信号以16Hz 抽样频率播放。 （2）用wavread 函数读取文件，对音乐信号进行内插或者抽取，再用sound 函数播放。

【仿真结果】

（1）32kHz 以16kHz 播放声音低沉，16kHz 能够正常播放，8kHz 声音变得尖细。 （2）音乐信号均正常播放。

【结果分析】

32kHz 音乐信号2倍抽取后，时域抽样间隔变为2倍，抽样频率变为原来二分之一，可以以16kHz

正常播放；8kHz音乐信号2倍内插，时域抽样间隔变为二分之一，抽样频率变为原来2倍，可以以16kHz正常播放。

【自主学习内容】

wavread函数与sound函数，内插函数及抽取函数。

【阅读文献】

【发现问题】(专题研讨或相关知识点学习中发现的问题)：

【问题探究】

【仿真程序】

（1）

x=wavread('C:\Users\Administrator\Desktop\数字信号研讨\kdqg32k.wav')

sound(x,16000);

x=wavread('C:\Users\Administrator\Desktop\数字信号研讨\kdqg16k.wav')

sound(x,16000);

x=wavread('C:\Users\Administrator\Desktop\数字信号研讨\kdqg8k.wav')

sound(x,16000);

（2）

x=wavread('C:\Users\Administrator\Desktop\数字信号研讨\kdqg32k.wav');

y=decimate(x,2);

sound(x,16000);

x=wavread('C:\Users\Administrator\Desktop\数字信号研讨\kdqg8k.wav');

y=interp(x,2);

sound(x,16000);

【研讨题目】实践题

3．对于人体脉搏信号可以综合地反映人体的心脏器官和心血管系统的生理变化。对人体脉搏信号进行非线性分析的研究，为评估人体心脏器官与血液循环系统的生理状态以及研究人的情绪对生理变化的影响提供了依据，并且对进一步深入认识人体生命系统的复杂性现象及规律具有重要意义。本题目将给出实际采集的人体脉搏信号，采样率为f sam=200Hz，

疲劳状态下脉搏信号：疲劳.txt

清晰状态下脉搏信号：清醒.txt

(1) 利用Matlab读取脉搏信号，并进行脉搏信号的频谱分析，得出人体脉搏频谱特征。

(2) 已知人体脉搏信号通常分布在1-20Hz内，设计合适的滤波器得到纯净的脉搏信号。

(3) 由于脉搏信号通常分布在1-20Hz内，而人耳听觉范围在300-3400Hz范围，为使人耳能够听到清醒的脉搏信号，请设计合适的多速率信号处理系统实现。

【频谱特征分析】

频谱主要部分是1-20Hz的低频分量。

【设计步骤】

1、先画出脉搏信号的时域波形和频谱。

2、通过级联的巴特沃斯滤波器滤波

3、通过多速率信号处理系统实现对脉搏信号的播放

【仿真结果】

（1）

（2）

（3）

【结果分析】

1、采用的脉搏信号中主要是低于20Hz的低频分量，在滤波是采用巴特沃斯高通滤波器和巴特沃斯低通滤波器的级联，比带通滤波器实现的代价要低。

2、通过内插滤波器可以实现对低频率信号的播放，内插产生镜像分量，取出可听的频率分量。

【仿真程序】

（1）

fs=200;

T=1/fs;

x1=load('C:\Users\Administrator\Desktop\数字信号研讨\清醒.txt');

x2=load('C:\Users\Administrator\Desktop\数字信号研讨\疲劳.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('频率');ylabel('幅度');

axis([0,100,0,3.5*10^5]);

subplot(2,2,4);plot(W2,abs(X2));

title('疲劳状态脉搏信号频谱');xlabel('频率');ylabel('幅度');

axis([0,100,0,3.5*10^5]);

（2）

fs=200;T=1/fs;

x=load('C:\Users\Administrator\Desktop\数字信号研讨\清醒.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('频率');ylabel('幅度'); axis([0,100,0,15000]);

fs=200;T=1/fs;

x=load('C:\Users\Administrator\Desktop\数字信号研讨\疲劳.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('频率');ylabel('幅度'); axis([0,100,0,15000]);

（3）

%以下通过内插实现对脉搏信号的播放

%先对人清醒状态的脉搏信号进行带通滤波，滤掉噪声信号

fs=200;T=1/fs;

x=load('C:\Users\Administrator\Desktop\数字信号研讨\清醒.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低通滤波器

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('幅度');

%构件带通滤波器

Ws1=0.01;Wp1=0.03;Wp2=0.1;Ws2=0.12;

ds=0.0017;dp=ds;

W1=[Ws1 Wp1 Wp2 Ws2];a=[0 1 0];dev=[ds dp ds]; [M,f0,a0,W0]=remezord(W1,a,dev);

h=remez(M,f0,a0,W0);

%滤除镜像部分

y0=conv(y1,h);

Y0=freqz(y0,1,W);

figure(3);plot(W/pi,abs(Y0));

title('滤除其他镜像部分后脉搏信号的频谱');

xlabel('归一化频率');ylabel('幅度');

sound(y0);