现代信号处理.docx
《现代信号处理.docx》由会员分享，可在线阅读，更多相关《现代信号处理.docx（21页珍藏版）》请在冰豆网上搜索。
现代信号处理
现代信号处理作业
班级:
姓名：
流水号：
1、设采样周期T=250μs（采样频率fs=4kHz），用脉冲响应不变法和双线性变换法设计一个三阶巴特沃兹滤波器，其3dB边界频率为fc=1kHz。
答：
MATLAB程序如下：
[B,A]=butter（3,2*pi*1000,'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）,'-'）;
grid;
title（‘现代信号处理第一题三阶巴特沃兹滤波器’）;
xlabel（'频率/Hz’）;
ylabel（'幅值/dB'）;
程序运行结果如下：
程序中第一个butter的边界频率2π×1000，为脉冲响应不变法原型低通滤波器的边界频率；第二个butter的边界频率2/T=2/0.00025，为双线性变换法原型低通滤波器的边界频率.图1给出了这两种设计方法所得到的频响，虚线为脉冲响应不变法的结果；实线为双线性变换法的结果。
脉冲响应不变法由于混叠效应，使得过渡带和阻带的衰减特性变差，并且不存在传输零点。
同时，也看到双线性变换法，在z=-1即Ω=π或f=2000Hz处有一个三阶传输零点，这个三阶零点正是模拟滤波器在ω=∞处的三阶传输零点通过映射形成的。
2、设计一数字高通滤波器，它的通带为400～500Hz，通带内容许有0.5dB的波动，阻带内衰减在小于317Hz的频带内至少为19dB，采样频率为1000Hz。
答：
MATLAB程序如下：
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;
title（‘现代信号处理第二题数字高通滤波器’）;
xlabel（'频率/Hz'）;
ylabel（'幅度/dB'）;
程序运行结果如下：
3、设计一巴特沃兹带通滤波器，其3dB边界频率分别为f2=110kHz和f1=90kHz，在阻带f3=120kHz处的最小衰减大于10dB，采样频率fs=400kHz。
答：
MATLAB程序如下：
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],[0.01wr],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;
title（‘现代信号处理第三题巴特沃兹带通滤波器’）
xlabel（'频率/kHz'）
ylabel（'幅度/dB'）
程序运行结果如下：
4、一数字滤波器采样频率fs=1kHz，要求滤除100Hz的干扰，其3dB的边界频率为95Hz和105Hz。
设计原型归一化滤波器。
答：
MATLAB程序如下：
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;
title（‘现代信号处理第四题归一化滤波器’）;
xlabel（'频率/Hz'）;
ylabel（'幅度/dB'）;
程序运行结果如下：
5、用计算机麦克录自己的语音信号，语音信号采样频率为22050，用matlab完成下列分析：
1）播放语音信号；对信号做1024点FFT变换；做原始语音信号的时域图形；绘制原始语音信号的频率响应图；做原始语音信号的FFT频谱图。
2）在语音信号中加入随机噪声；播放加噪声后的语音信号，绘制加噪后的语音信号；
3）设计合适的数字滤波器，将上述加噪声滤掉；播放滤波后的信号；绘制滤波前和滤波后的语音信号及频谱图。
答：
1）MATLAB程序如下：
fs=22050;
[x1,fs,bits]=wavread（'D:
\11.wav'）;
sound（x1,fs,bits）;
figure
（1）;
plot（x1）;
title（'原始信号波形图'）;
xlabel（'时间轴'）;
ylabel（'幅值'）;
figure
（2）;
freqz（x1）;
title（'原始语音信号频率响应图'）;
y1=fft（x1,1024）;
f=fs*（0:
511）/1024;
figure（3）;
plot（abs（y1（1:
512）））%做原始语音信号的FFT频谱图
title（'原始语音信号FFT频谱'）;
1）程序运行结果如下：
2）MATLAB程序如下：
fs=22050;
[y,fs,nbits]=wavread（'D:
\11.wav'）;
n=length（y）;%求出语音信号的长度
Noise=0.01*randn（n,2）;%随机函数产生噪声
Si=y+Noise;%语音信号加入噪声
sound（Si,fs,nbits）;%回放语音信号
subplot（2,1,1）;
plot（Si）;title（'加噪语音信号的时域波形'）;
y1=fft（y,1024）;
f=fs*（0:
511）/1024;
subplot（2,1,2）;
plot（abs（y（1:
512）））%做原始语音信号的FFT频谱图
title（'加噪语音信号的FFT频谱'）;
2）程序运行结果如下：
3）MATLAB程序如下：
fs=22050;
[y,fs,nbits]=wavread（'D:
\11.wav'）;
n=length（y）;%求出语音信号的长度
Noise=0.01*randn（n,2）;%随机函数产生噪声
x2=y+Noise;%语音信号加入噪声
wp=0.25*pi;
ws=0.3*pi;
Rp=1;
Rs=15;
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）;%绘制频率响应曲线
plot（W*fs/（2*pi）,abs（H））
grid
xlabel（'频率／Hz'）
ylabel（'频率响应幅度'）
title（'Butterworth'）
f1=filter（bz,az,x2）;
figure
（1）
subplot（2,1,1）
plot（x2）%画出滤波前的时域图
title（'滤波前的时域波形'）;
xlabel（'频率/Hz'）;
ylabel（'幅值'）;
subplot（2,1,2）
plot（f1）;%画出滤波后的时域图
title（'滤波后的时域波形'）;
xlabel（'频率/Hz'）;
ylabel（'幅值'）;
sound（f1,22050）;%播放滤波后的信号
F0=fft（f1,1024）;
f=fs*（0:
511）/1024;
figure
（2）
y2=fft（x2,1024）;
subplot（2,1,1）;
plot（f,abs（y2（1:
512）））;%画出滤波前的频谱图
title（'滤波前的频谱'）
xlabel（'频率/Hz'）;
ylabel（'幅值'）;
subplot（2,1,2）
F1=plot（f,abs（F0（1:
512）））;%画出滤波后的频谱图
title（'滤波后的频谱'）
xlabel（'频率/Hz'）;
ylabel（'幅值'）;
3）程序运行结果如下：
6、信号
，信号带宽
，采样频率
，采样点数
。
用STFT和WD分析其特性。
答：
1）STFT分析，MATLAB程序如下：
k=4;T=6;
fc=k*T;
fs=4*fc;
Ts=1/fs;
N=T/Ts;
x=zeros（1,N）;
t=0:
N-1;
x=exp（j*k*pi*（t*Ts）.^2）;
subplot（2,2,1）;
plot（t*Ts,real（x））;
title（'原信号'）;
X=fft（x）;
X=fftshift（X）;
subplot（2,2,2）;
plot（（t-N/2）*fs/N,abs（X））;
title（'FFT的结果'）;
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））;
title（'STFT的结果'）;
subplot（2,2,4）;
contour（F,T,abs（TF））;
title（'STFT的俯视图'）;
1）STFT分析，程序运行结果如下：
2）WD分析，MATLAB程序如下：
k=4;T=6;
fc=k*T;
fs=4*fc;
Ts=1/fs;
N=T/Ts;
x=zeros（1,N）;
t=0:
N-1;
x=exp（j*k*pi*（t*Ts）.^2）;
subplot（2,2,1）;
plot（t*Ts,real（x））;
title（'原信号'）;
X=fft（x）;
X=fftshift（X）;
subplot（2,2,2）;
plot（（t-N/2）*fs/N,abs（X））;
title（'FFT的结果'）;
R=zeros（N,N）;
forn=0:
N-1,
M=min（n,N-1-n）;
fork=0:
M,
R（n+1,k+1）=x（n+k+1）*conj（x（n-k+1））;
end
fork=N-1:
-1:
N-M,
R（n+1,k+1）=conj（R（n+1,N-k+1））;
end
end
TF=zeros（N,N）;
forn=0:
N-1,
temp=fft（R（n+1,:
））;
temp=fftshift（temp）;
TF（n+1,:
）=temp;
end
fnew=（t-N/2）*fs/2/N;
tnew=（0:
N-1）*Ts;
[F,T]=meshgrid（fnew,tnew）;
subplot（2,2,3）;
mesh（F,T,abs（TF））;
title（'WD的结果'）;
subplot（2,2,4）
contour（F,T,abs（TF））;
title（'WD的俯视图'）;
2）WD分析，程序运行结果如下：
7、数字信号X（n）,0≤n≤N-1,N=1000,在（125，250）和（500，625）分别有两个频率不同的正弦信号，用STFT和WD分析其特性。
答：
1）STFT分析，MATLAB程序如下：
N=1000;t=0:
N-1;
x=zeros（1,N）;
x（125:
250）=cos（pi*（t（125:
250）-125）/10）;
x（500:
625）=cos（pi*（t（500:
625）-500）/5）;
subplot（2,2,1）;
plot（x）;
title（'原信号'）;
X=fft（x）;
X=fftshift（X）;
subplot（2,2,2）;
plot（abs（X））;
title（'FFT的结果'）;
Nw=20;
L=Nw/2;
Ts=（N-Nw）/L+1;
nfft=32;
TF=zeros（Ts,nfft）;
fori=1:
Ts,
xw=x（（i-1）*L+1:
i*L+L）;
temp=fft（xw,nfft）;
temp=fftshift（temp）;
TF（i,:
）=temp;
end
subplot（2,2,3）;
mesh（abs（TF））;
title（'STFT的结果'）;
subplot（2,2,4）;
contour（abs（TF））;
title（'STFT的俯视图'）;
1）STFT分析，程序运行结果如下：
2）WD分析，MATLAB程序如下：
N=1000;t=0:
N-1;
x=zeros（1,N）;
x（125:
250）=cos（pi*（t（125:
250）-125）/10）;
x（500:
625）=cos（pi*（t（500:
625）-500）/5）;
subplot（2,2,1）;
plot（t,real（x））;
title（'原信号'）;
X=fft（x）;
X=fftshift（X）;
subplot（2,2,2）;
plot（（t-N/2）/N,abs（X））;
title（'FFT的结果'）;
R=zeros（N,N）;
forn=0:
N-1,
M=min（n,N-1-n）;
fork=0:
M,
R（n+1,k+1）=x（n+k+1）*conj（x（n-k+1））;
end
fork=N-1:
-1:
N-M,
R（n+1,k+1）=conj（R（n+1,N-k+1））;
end
end
TF=zeros（N,N）;
forn=0:
N-1,
temp=fft（R（n+1,:
））;
temp=fftshift（temp）;
TF（n+1,:
）=temp;
end
fnew=（t-N/2）/N;
tnew=0:
N-1;
[F,T]=meshgrid（fnew,tnew）;
subplot（2,2,3）;
mesh（F,T,abs（TF））;
title（'WD的结果'）;
subplot（2,2,4）
contour（F,T,abs（TF））;
title（'WD的俯视图'）;
2）WD分析，程序运行结果如下：
升级会员 