信号处理作业.docx
《信号处理作业.docx》由会员分享，可在线阅读，更多相关《信号处理作业.docx（15页珍藏版）》请在冰豆网上搜索。
信号处理作业
现代信号处理
学院：
仪器与电子学院
专业：
仪器仪表工程
班级：
姓名：
学号：
流水号：
时间：
2014年6月27日
1、设采样周期T=250μs（采样频率fs=4kHz），用脉冲响应不变法和双线性变换法设计一个三阶巴特沃兹滤波器，其3dB边界频率为fc=1kHz。
设计步骤如下：
（1）确定所需类型数字滤波器的技术指标。
（2）将所需类型数字滤波器的边界频率转换成相应的模拟滤波器的边界频率，转换公式为Ω=2/Ttan（0.5ω）
（3）将相应类型的模拟滤波器技术指标转换成模拟低通滤波器的技术指标。
（4）设计模拟低通滤波器。
（5）通过频率变换将模拟低通转换成相应类型的过渡模拟滤波器。
（6）采用脉冲响应不变法和双线性变换法将相应类型的过渡模拟滤波器转换成所需类型的数字滤波器。
程序设计如下：
[B,A]=butter（3,2*pi*1000,'s'）;%第一个参数为阶数；第二个参数为角频率；第三个参数为系数B、A按s的正降幂排列。
[num1,den1]=impinvar（B,A,4000）;%用冲激响应不变法将模拟滤波器变换成数字滤波器
[h1,w]=freqz（num1,den1）;%获取滤波器的频率特性
[B,A]=butter（3,2/0.00025,'s'）;
[num2,den2]=bilinear（B,A,4000）;%实现双线性变换，即由模拟滤波器得到数字滤波器
[h2,w]=freqz（num2,den2）;
f=w/pi*2000;
plot（f,abs（h1）,'-.',f,abs（h2）,'-'）;% x为图形上之x坐标向量，y为其对应的y坐标向量
grid;
xlabel（'频率/Hz'）
ylabel（'幅值/dB'）
程序中第一个butter的边界频率2π×1000，为脉冲响应不变法原型低通滤波器的边界频率；第二个butter的边界频率2/T=2/0.00025，为双线性变换法原型低通滤波器的边界频率.图1给出了这两种设计方法所得到的频响，虚线为脉冲响应不变法的结果；实线为双线性变换法的结果。
脉冲响应不变法由于混叠效应，使得过渡带和阻带的衰减特性变差，并且不存在传输零点。
同时，也看到双线性变换法，在z=-1即Ω=π或f=2000Hz处有一个三阶传输零点，这个三阶零
点正是模拟滤波器在ω=∞处的三阶传输零点通过映射形成的。
MATLAB计算结果如图1所示：
图1
2、设计一数字高通滤波器，它的通带为400～500Hz，通带内容许有0.5dB的波动，阻带内衰减在小于317Hz的频带内至少为19dB，采样频率为1000Hz。
wc=2*1000*tan（2*pi*400/（2*1000））;
wt=2*1000*tan（2*pi*317/（2*1000））;
[N,wn]=cheb1ord（wc,wt,0.5,19,'s'）;
[B,A]=cheby1（N,0.5,wn,'high','s'）;
[num,den]=bilinear（B,A,1000）;
[h,w]=freqz（num,den）;
f=w/pi*500;
plot（f,20*log10（abs（h）））;
axis（[0,500,-80,10]）;
grid;
xlabel（''）
ylabel（'幅度/dB'）
图2给出了MATLAB计算的结果，可以看到模拟滤波器在Ω=∞处的三阶零点通过高通变换后出现在ω=0（z=1）处，这正是高通滤波器所希望得到的。
图2
3、设计一巴特沃兹带通滤波器，其3dB边界频率分别为f2=110kHz和f1=90kHz，在阻带f3=120kHz处的最小衰减大于10dB，采样频率fs=400kHz。
w1=2*400*tan（2*pi*90/（2*400））;
w2=2*400*tan（2*pi*110/（2*400））;
wr=2*400*tan（2*pi*120/（2*400））;
[N,wn]=buttord（[w1w2],[1wr],3,10,'s'）;
[B,A]=butter（N,wn,'s'）;
[num,den]=bilinear（B,A,400）;
[h,w]=freqz（num,den）;
f=w/pi*200;
plot（f,20*log10（abs（h）））;
axis（[40,160,-30,10]）;
grid;
xlabel（'频率/kHz'）
ylabel（'幅度/dB'）
MATLAB计算结果如图3所示：
图3
4、一数字滤波器采样频率fs=1kHz，要求滤除100Hz的干扰，其3dB的边界频率为95Hz和105Hz。
低通滤波器传递函数如下：
H1a（s）=1/（1+s）
w1=95/500;
w2=105/500;
[B,A]=butter（1,[w1,w2],'stop'）;
[h,w]=freqz（B,A）;
f=w/pi*500;
plot（f,20*log10（abs（h）））;
axis（[50,150,-30,10]）;
grid;
xlabel（'频率/Hz'）
ylabel（'幅度/dB'）
图4为MATLAB的计算结果
图4
5、用计算机麦克录自己的语音信号，语音信号采样频率为22050，用matlab完成下列分析：
1）播放语音信号；对信号做1024点FFT变换；做原始语音信号的时域图形；绘制原始语音信号的频率响应图；做原始语音信号的FFT频谱图。
fs=22050;%语音信号采样频率为22050
x1=wavread（'D:
\program-tool\MATLAB\R2010b\work\clock001.wav'）;%读取语音信号的数据，赋给变量x1
sound（x1,32768）;%播放语音信号
y1=fft（x1,1024）; %对信号做1024点FFT变换
f=fs*（0:
511）/1024;
figure
（1）;
plot（x1）; %做原始语音信号的时域图形
title（'原始语音信号'）;
xlabel（'timen'）;
ylabel（'fuzhin'）;
figure
（2）;
freqz（x1）;%绘制原始语音信号的频率响应图
title（'频率响应图'）;
figure（3）;
subplot（2,1,1）;
plot（abs（y1（1:
512）））;%做原始语音信号的FFT频谱图
title（'原始语音信号FFT频谱'）;
subplot（2,1,2）;
plot（f,abs（y1（1:
512）））;
title（'原始语音信号频谱'）;
xlabel（'Hz'）;
ylabel（'fuzhi'）;
MATLAB计算结果如图5、6、7所示：
图5
图6
图7
2）在语音信号中加入随机噪声；播放加噪声后的语音信号，绘制加噪后的语音信号；
首先画出语音信号的时域波形，再对语音信号进行快速傅里叶变换，得到信号的频谱特性。
给原始的语音信号加上一个高频余弦噪声，频率为5kHz。
画出加噪后的语音信号时域和频谱图，与原始信号对比，可以很明显的看出区别
程序如下：
fs=22050;
x1=wavread（'qq.wav'）;
f=fs*（0:
511）/1024;
t=0:
1/22050:
（size（x1）-1）/22050;%将所加噪声信号的点数调整到与原始信号相同
Au=0.3;
d=[Au*cos（2*pi*5000*t）]';%噪声为5kHz的余弦信号
x2=x1+d;
sound（x2,22050）;%播放加噪声后的语音信号
y2=fft（x2,1024）;
figure
（1）
plot（t,x2）;
title（'加噪后的信号'）;
xlabel（'timen'）;
ylabel（'fuzhin'）;
figure
（2）
subplot（2,1,1）;
plot（f,abs（y1（1:
512）））;
title（'原始语音信号频谱'）;
xlabel（'Hz'）;
ylabel（'fuzhi'）;
subplot（2,1,2）;
plot（f,abs（y2（1:
512）））;
title（'加噪后的信号频谱'）;
xlabel（'Hz'）;
ylabel（'fuzhi'）;
得到原始语音信号时域图形，原始语音信号幅度谱，加入噪声之后的语音信号时域图形，加入噪声之后的语音信号幅度谱如图8、9所示：
图8
图9
3）设计合适的数字滤波器，将上述加噪声滤掉；播放滤波后的信号；绘制滤波前和滤波后的语音信号及频谱图。
双线性变换法设计Butterworth滤波器
fs=22050;
x1=wavread（'qq.wav'）;
t=0:
1/22050:
（size（x1）-1）/22050;
Au=0.03;
d=[Au*cos（2*pi*10000*t）]';
x2=x1+d;
wp=0.25*pi;
ws=0.3*pi;
Rp=1;Rs=15;Fs=32768;Ts=1/Fs;
wp1=2/Ts*tan（wp/2）;%将模拟指标转换成数字指标
ws1=2/Ts*tan（ws/2）;
[N,Wn]=buttord（wp1,ws1,Rp,Rs,'s'）;%选择滤波器的最小阶数
[Z,P,K]=buttap（N）;%创建butterworth模拟滤波器
[Bap,Aap]=zp2tf（Z,P,K）;
[b,a]=lp2lp（Bap,Aap,Wn）;
[bz,az]=bilinear（b,a,Fs）;%用双线性变换法实现模拟滤波器到数字滤波器的转换
[H,W]=freqz（bz,az）;%绘制频率响应曲线
figure
（1）
plot（W*Fs/（2*pi）,abs（H））
grid
xlabel（'频率／Hz'）
ylabel（'频率响应幅度'）
title（'Butterworth'）
f1=filter（bz,az,x2）;
figure
（2）
subplot（2,1,1）
plot（t,x2）%画出滤波前的时域图
title（'滤波前的时域波形'）;
subplot（2,1,2）
plot（t,f1）;%画出滤波后的时域图
title（'滤波后的时域波形'）;
sound（f1,22050）;%播放滤波后的信号
F0=fft（f1,1024）;
f=fs*（0:
511）/1024;
figure（3）
y2=fft（x2,1024）;
subplot（2,1,1）;
plot（f,abs（y2（1:
512）））;%画出滤波前的频谱图
title（'滤波前的频谱'）
xlabel（'Hz'）;
ylabel（'fuzhi'）;
subplot（2,1,2）
F1=plot（f,abs（F0（1:
512）））;%画出滤波后的频谱图
title（'滤波后的频谱'）
xlabel（'Hz'）;
ylabel（'fuzhi'）;
计算结果如图10、11、12
图10
图11
图12
6、信号
，信号带宽
，采样频率
，采样点数
。
用STFT和WD分析其特性。
设计过程：
用短时傅里叶变换分析其在时域、频域特性，经过STFT和WD运算后的时频分布结果特点，从而将它们识别。
设计的源程序如下：
k=4;t=5;fc=k*T=20;fs=4*fc;Ts=1/fs;N=T/Ts=400;
X=zeros（1,N）;
T=0:
N-1;
X=exp（j*k*pi*（t*Ts）.＾2）;
Subplot（2,2,2）;
Plot（（t-N/2）*fs/N,abs（X））;
Nw=20;
L=Nw/2;
Tn=（N-Nw）/L+1;
Nfft=32;
TF=zeros（Tn,nfft）;
fori=1:
Tn
xw=x（（i-1）*10+1:
i*10+10）;
temp=fft（xw,nfft）;
temp=fftshift（temp）;
TF（i,:
）=temp;
end
subplot（2,2,3）;
fnew=（（1:
nfft）-nfft/2）*fs/nfft;
tnew=（1:
Tn）*L*Ts;
[F,T]=meshgrid（fnew,tnew）;
mesh（F,T,abs（TF）;
subplot（2,2,4）;
计算结果如图13所示
图13
7、数字信号
在（125，250）和（500，625）分别有两个频率不同的正弦信号，用STFT和WD分析其特性。
k=4;T=5;fc=20;fs=4*fc;Ts=1/fs;N=T/Ts;
x=zeros（1,N）;
t=125:
250;
x=sin（t*Ts）;
subplot（2,2,1）
plot（t*Ts,real（x））;
X=fft（x）;
X=fftshift（X）;
subplot（2,2,2）;
plot（（t-N/2）*fs/N,abs（X））;
t=500:
625;
x=sin（t*Ts）;
subplot（2,2,1）
plot（t*Ts,real（x））;
X=fft（x）;
X=fftshift（X）;
subplot（2,2,2）;
plot（（t-N/2）*fs/N,abs（X））;
计算结果如图14所示
图14
升级会员 