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

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

指导教师 薛健
时间 2014.6
【目的】
(1) 掌握IIR和FIR数字滤波器的设计和应用；
(2) 掌握多速率信号处理中的基本概念和方法 ；
(3) 学会用Matlab计算小波分解和重建。
(4)了解小波压缩和去噪的基本原理和方法。

【研讨题目】
一、
(1)播放音频信号 yourn.wav，确定信号的抽样频率，计算信号的频谱，确定噪声信号的频率范围；
(2)设计IIR数字滤波器，滤除音频信号中的噪声。通过实验研究，的选择对滤波效果及滤波器阶数的影响，给出滤波器指标选择的基本原则，确定你认为最合适的滤波器指标。
（3）设计FIR数字滤波器，滤除音频信号中的噪声。与（2）中的IIR数字滤波器，从滤波效果、幅度响应、相位响应、滤波器阶数等方面进行比较。
【设计步骤】

【仿真结果】

【结果分析】
由频谱知噪声频率大于3800Hz。FIR和IIR都可以实现滤波，但从听觉上讲，人对于听觉不如对图像（视觉）明感，没必要要求线性相位，因此，综合来看选IIR滤波器好一点，因为在同等要求下，IIR滤波器阶数可以做的很低而FIR滤波器阶数太高，自身线性相位的良好特性在此处用处不大。
【自主学习内容】
MATLAB滤波器设计
【阅读文献】
老师课件，教材
【发现问题】 (专题研讨或相关知识点学习中发现的问题)：
过渡带的宽度会影响滤波器阶数N
【问题探究】
通过实验，但过渡带越宽时，N越小，滤波器阶数越低，过渡带越窄反之。这与理论相符合。
【仿真程序】
信号初步处理部分：
[x1,Fs,bits] = wavread('yourn.wav');
sound(x1,Fs);
y1=fft(x1,1024);
f=Fs*(0:511)/1024;
figure(1)
plot(x1)
title('原始语音信号时域图谱');
xlabel('time n');
ylabel('magnitude n');
figure(2)
freqz(x1)
title('频率响应图')
figure(3)
subplot(2,1,1);
plot(abs(y1(1:512)))
title('原始语音信号FFT频谱')
subplot(2,1,2);
plot(f,abs(y1(1:512)));
title(‘原始语音信号频谱')
xlabel('Hz');
ylabel('magnitude');
IIR：
fp=2500;fs=3500;
wp = 2*pi*fp/FS;
ws = 2*pi*fs/FS;
Rp=1;
Rs=15;
Ts=1/Fs;
wp = 2*pi*fp/FS;
ws = 2*pi*fs/FS;
wp1=2/Ts*tan(wp/2);
ws1=2/Ts*tan(ws/2);
t=0:1/11000:(size(x1)-1)/11000;
Au=0.03;
d=[Au*cos(2*pi*5000*t)]';
x2=x1+d;

[N,Wn]=buttord(wp1,ws1,Rp,Rs,'s');
[Z,P,K]=buttap(N);
[Bap,Aap]=zp2tf(Z,P,K);
[b,a]=lp2lp(Bap,Aap,Wn);
[bz,az]=bilinear(b,a,Fs); %
[H,W]=freqz(bz,az);
figure(4)
plot(W*Fs/(2*pi),abs(H))
grid
xlabel('频率/Hz')
ylabel('频率响应幅度')
title('Butterworth')
f1=filter(bz,az,x2);
figure(5)
subplot(2,1,1)
plot(t,x2)
title('滤波前时域波形');
subplot(2,1,2)
plot(t,f1);
title('滤波后时域波形');
sound(f1,FS);
FIR
[x1,Fs,bits] = wavread('I:/dsp_2014_project3/yourn');
fp=2500;fs=3500;
wp = 2*pi*fp/Fs;
ws=2*pi*fs/Fs; Rs=50;
M=ceil((Rs-7.95)/(ws-wp)/2.285);
M=M+mod(M,2);
beta=0.1102*(Rs-8.7);
w=kaiser(M+1,beta);
wc=(wp+ws)/2;
alpha=M/2;
k=0:M;
hd=(wc/pi)*sinc((wc/pi)*(k-alpha));
h=hd.*w';
f1=filter(h,[1],x1);
[mag,W]=freqz(h,[1]);
figure(1)
plot(W*Fs/(2*pi),abs(mag));
grid;
xlabel('频率/Hz');
ylabel('频率响应幅度');
title('Kaiser´窗设计Ⅰ型线性相位FIR低通滤波器');
figure(2)
subplot(2,1,1)
plot(t,x1)
title('滤波前时域波形');
subplot(2,1,2)
plot(t,f1);
title('滤波后时域波形');
sound(f1,Fs);

二、(1)音频信号kdqg24k.wav抽样频率为24kHz，用
y = wavread('kdqg24k');
sound(y,16000);
播放该信号。试用频域的方法解释实验中遇到的现象;
(2)设计一数字系统 ，使得sound(y,16000)可播放出正常的音频信号；讨论滤波器的频率指标、滤波器的的类型(IIR,FIR)对系统的影响。
【仿真结果】

【结果分析】
24K的信号用16K播放，频谱会被拉宽，无法正常播放，通过2倍内插，通过滤波器，然后3倍抽取，得到的信号用16K播放器就能正常播放。
【自主学习内容】
功能：对时间序列进行重采样。
格式：
1.y = resample(x, p, q)
采用多相滤波器对时间序列进行重采样，得到的序列y的长度为原来的序列x的长度的p/q倍，p和q都为正整数。此时，默认地采用使用FIR方法设计的抗混叠的低通滤波器。
2.y = resample(x, p, q, n)
采用chebyshevIIR型低通滤波器对时间序列进行重采样，滤波器的长度与n成比例，n缺省值为10.
3.y = resample(x, p, q, n, beta)
beta为设置低通滤波器时使用Kaiser窗的参数，缺省值为5.
4.y = resample(x, p, q, b)
b为重采样过程中滤波器的系数向量。
5.[y, b] = resample(x, p, q)
输出参数b为所使用的滤波器的系数向量。
说明：
x--时间序列
p、q--正整数，指定重采样的长度的倍数。
n--指定所采用的chebyshevIIR型低通滤波器的阶数，滤波器的长度与n成比列。
beta--设计低通滤波器时使用Kaiser窗的参数，缺省值为5.

【阅读文献】
PPt 课本
【发现问题】 (专题研讨或相关知识点学习中发现的问题)：
采样频率与播放频率之间不是整数倍关系
【问题探究】
此时内插和抽取结合实现正常播放
【仿真程序】
fs=24000;
x1= wavread('I:/dsp_2014_project3/kdqg24k.wav');
sound(x1,16000);
y1=fft(x1,1024);
f1=fs*(0:511)/1024;
f2=fs/2*3*(0:511)/1024;

figure(1)
subplot(2,1,1);
plot(f1,abs(y1(1:512)));
title('原始语音信号24K正常播放频谱')
xlabel('Hz');
ylabel('magnitude');

subplot(2,1,2);
plot(f2,abs(y1(1:512)));
title('原始语音信号16K播放频谱')
xlabel('Hz');
ylabel('magnitude');
y = resample(x1,2, 3);
sound(y,16000);
y2=fft(y,1024);
figure(2)
subplot(2,1,1);
plot(f1,abs(y1(1:512)));
title('原始语音信号24K正常播放频谱')
xlabel('Hz');
ylabel('magnitude');

subplot(2,1,2);
plot(f2,abs(y2(1:512)));
title('原始语音信号16K经过设计的数字系统后播放频谱')
xlabel('Hz');
ylabel('magnitude');
三、 对连续信号
x(t)=40t2(1-t)4cos(12pt)[0<t<1] +40(t-1)4(2-t)2cos(40pt)[1<t<2]+0.1n(t)
在区间[0，2]均匀抽样1024点得离散信号x[k]，其中n(t)是零均值方差为1的高斯噪声
(1)画出信号x(t)的波形；
(2)计算并画出db7小波的5级小波变换系数；
(3)通过观察小波系数，确定阈值化处理的阈值；
(4)对小波系数进行阈值化处理，画出去噪后的信号波形，求出最大误差和均方误差；
(5)对近似系数和小波系数均进行阈值化处理，画出去噪后的信号波形，求出最大误差和均方误差；
(6)用Haar小波基，重复(3)-(5)；
(7)讨论所得结果。
【仿真结果】
（1）

（2）

（3）T=4.3000 t=96.84%
(4)db7小波基

Emabs=0.9166
(5)Haar 小波基 T=3.100; t=96.18%

Emabs=1.1341
(6)db14小波基 T=4.6，t=94.39%

Emabs=0.8965

【结果分析】
选用dbp系列小波基对信号去噪时，p值越大，去噪的效果越好，最大重建误差也小
【自主学习内容】
MATLAB小波变换
【阅读文献】
PPT和课本
【发现问题】 (专题研讨或相关知识点学习中发现的问题)：
选用dbp系列小波基对信号去噪时，p值越大，去噪的效果越好，最大重建误差也小
【问题探究】
不同dp小波基不同，会对结果产生不同影响
【仿真程序】
N=4096;
k=linspace(0,2,N);
nt=randn(size(k));
x=40*k.^2.*(1-k).^4.*cos(12*pi*k).*(0<k&k<1)+40.*(k-1).^4.*(2-k).^2.*cos(80*pi*k).*(1<k&k<2)+0.1*nt;
figure;plot(x);title('signal with noise');

dwtmode('per');
[C,L] = wavedec(x,6,'db14')
figure;plot(k,C);title('Wavelet coefficients');

M=0;for k=1:4096;T=4.6;
if abs(C(1,k))<=T;C(1,k)=0;end
if C(1,k)~=0;M=M+1; end
end

A1=C.*C;U1=0;A2=x.*x;U2=0;
for k=1:1024;
U1=U1+A1(1,k);
U2=U2+A2(1,k);
end
t=U1/U2

[XD,CXD,LXD]=wden(x,'heursure','s','one',6,'db14');
s=waverec(CXD,LXD,'db14');
figure;subplot(211);plot(s);title('Ronconstructed signal with no niose');
subplot(212);
plot(x-s);title('Error');
d=max(x-s);