数字信号.docx
《数字信号.docx》由会员分享，可在线阅读，更多相关《数字信号.docx（20页珍藏版）》请在冰豆网上搜索。
数字信号
附件一：
%巴特沃斯模拟高通滤波器
%阻带截止频率fp=1.5kHz，阻带内最大衰减ap=2dB；通带截止频率fs=2kHz，通带最小衰减as=30dB。
wp=2*pi*1500;ws=2*pi*2000;Rp=2;As=30;
[N,wc]=buttord（wp,ws,Rp,As,'s'）;%计算滤波器阶数N和3dB截止频率wc
[B,A]=butter（N,wc,'high','s'）;%计算滤波器系统函数分子分母多项式系数
k=0:
511;fk=0:
14000/512:
14000;wk=2*pi*fk;
Hk=freqs（B,A,wk）;
figure
（2）;
plot（fk/1000,20*log10（abs（Hk）））;
xlabel（'频率（kHz）'）;ylabel（'幅度（dB）'）
title（'巴特沃斯模拟滤波器的频率响应'）;
axis（[0,14,-40,5]）
gridon
附件二
%椭圆模拟低通滤波器
%通带截止频率fp=1.5kHz，通带内最大衰减ap=2dB；阻带截止频率fs=2kHz，阻带最小衰减as=30dB。
wp=2*pi*1500;ws=2*pi*2000;Rp=2;As=30;
[N,wpo]=ellipord（wp,ws,Rp,As,'s'）;%计算滤波器阶数N和3dB截止频率wc
[B,A]=ellip（N,Rp,As,wpo,'s'）;%计算滤波器系统函数分子分母多项式系数
k=0:
511;fk=0:
14000/512:
14000;wk=2*pi*fk;
Hk=freqs（B,A,wk）;
figure
（1）;
plot（fk/1000,20*log10（abs（Hk）））;
xlabel（'频率（kHz）'）;ylabel（'幅度（dB）'）
title（'椭圆模拟滤波器的频率响应'）;
axis（[0,10,-40,5]）
gridon
附件三
%巴特沃斯数字高通滤波器
[x,fs,bits]=wavread（'gzw.wav'）;%播放原始信号
N=length（x）;%返回采样点数
t=（1:
N）/fs;
df=fs/N;%采样间隔
n1=1:
N/2;
f=（n1-1）*df;%频带宽度
%************************巴特沃斯滤波器设计*****************************
FS=8000;
%通带、阻带截止频率
Fl=1500;Fh=2000;
%频率预畸
wp=（Fl/FS）*2*pi;%临界频率采用角频率表示
ws=（Fh/FS）*2*pi;%临界频率采用角频率表示
OmegaP=2*FS*tan（wp/2）;
OmegaS=2*FS*tan（ws/2）;
[k,Wn]=buttord（OmegaP,OmegaS,0.1,60,'s'）;
[b,a]=butter（k,Wn,'high','s'）;
[bz,az]=bilinear（b,a,FS）;%双线形变换映射为数字的
%绘制结果
H=freqz（bz,az,FS,'whole'）;
figure
（1）;
plot（abs（H））;
title（'巴特沃斯滤波器的频率响应'）;
xlabel（'频率/Hz'）;
ylabel（'幅值'）;
%**********************************************************************
ys=filter（bz,az,x）;%信号送入滤波器滤波，ys为输出
fftwave=fft（ys）;%将滤波后的语音信号进行快速傅立叶变换
figure
（2）;
plot（f,abs（fftwave（n1）））;
title（'巴特沃斯滤波器滤波后信号的频谱图'）;
xlabel（'频率/Hz'）;
ylabel（'幅值/V'）;
grid;
附件四
%椭圆数字低通滤波器
[x,fs,bits]=wavread（'gzw.wav'）;%播放原始信号
N=length（x）;%返回采样点数
t=（1:
N）/fs;
df=fs/N;%采样间隔
n1=1:
N/2;
f=（n1-1）*df;%频带宽度
FS=8000;
%通带、阻带截止频率
Fl=1500;Fh=2000;
%频率预畸
wp=（Fl/FS）*2*pi;%临界频率采用角频率表示
ws=（Fh/FS）*2*pi;%临界频率采用角频率表示
OmegaP=2*FS*tan（wp/2）;
OmegaS=2*FS*tan（ws/2）;
[k,Wn]=ellipord（OmegaP,OmegaS,0.2,60,'s'）;
[b,a]=ellip（k,0.2,60,Wn,'s'）;
[bz,az]=bilinear（b,a,FS）;%双线形变换映射为数字的
%绘制结果
H=freqz（bz,az,FS,'whole'）;
figure
（1）;
plot（abs（H））;
title（'椭圆滤波器的频率响应'）;
xlabel（'频率/Hz'）;
ylabel（'幅值'）;
%**********************************************************************
ys=filter（bz,az,x）;%信号送入滤波器滤波，ys为输出
fftwave=fft（ys）;%将滤波后的语音信号进行快速傅立叶变换
figure
（2）;
plot（f,fftwave（n1））;
title（'椭圆滤波器滤波后信号的频谱图'）;
xlabel（'频率/Hz'）;
ylabel（'幅值/V'）;
grid;
附件五：
clearall;
[x,Fs,bits]=wavread（'E:
/gzw/gzw.wav'）;%播放原始信号
N=length（x）;%返回采样点数
t=（1:
N）/Fs;
df=Fs/N;%采样间隔
n1=1:
N/2;
f=（n1-1）*df;%频带宽度
figure
（1）;
%将语音数据通过图像显示出来
plot（x,'LineWidth',2）,gridon;%信号的时域波形
title（'原始信号的时域波形'）;
xlabel（'时间/t'）;
ylabel（'幅值/A'）;
y0=fft（x）;%快速傅立叶变换
figure
（2）;
plot（f,20*log10（abs（y0（n1））））,gridon;%离散信号的频谱图
title（'原始信号的频谱图'）;
xlabel（'频率/Hz'）;
ylabel（'幅值/db'）;
grid;
%矩形窗设计低通滤波器，
%通带截止频率1.5khz,阻带截止频率2khz,阻带最小衰减20
fp=1500;fs=2000;rs=20;
%转换为数字指标
wp=2*pi*fp/Fs;
ws=2*pi*fs/Fs;
%计算频带宽度
Bt=ws-wp;
wc=（wp+ws）/2/pi;
%计算窗长度
M=ceil（1.8*pi/Bt）;
%利用fir1进行矩形滤波器设计
hn=fir1（M,wc,boxcar（M+1））;
figure（3）;
%显示矩形低通滤波器
stem（boxcar（M+1））;
title（'矩形低通滤波器'）;
figure（4）;
%计算滤波器的频响
[H,f]=freqz（hn,1,512,Fs）;
mag=20*log10（abs（H））;
plot（f,mag）;
title（'矩形低通频率响应'）;
%利用矩形窗进行滤波
yout=fftfilt（hn,x）;
figure（5）;
%显示滤波后的信号
plot（yout）;
title（'矩形低通滤波后的语音信号'）;
figure（6）;
%滤波后信号FFT变换
FYLB=fft（yout）;
FABS=abs（FYLB（1:
N/2））;
%显示滤波后信号的频谱
plot（FABS）;
title（'矩形低通滤波后的频谱'）;
%用汉宁窗进行高通滤波通带最大衰减1db，阻带最小衰减30db
%通带截止频率2.5khz
wp=0.5*pi;
%阻带截止频率1.5khz
ws=0.375*pi;
%计算频带宽度
B=wp-ws;
%计算窗函数宽度
N_h=ceil（6.2*pi/B）;
N_hanning=N_h+mod（N_h+1,2）;
wc=（wp+ws）/2/pi;
%调用fir1设计hanning高通数字滤波器
hn=fir1（N_hanning-1,wc,'low',hanning（N_hanning））;
figure（7）;
%显示汉宁滤波器
title（'汉宁滤波'）;
stem（hanning（N_hanning））;
figure（8）;
%计算汉宁滤波器的频响
[H,f]=freqz（hn,1,512,Fs）;
mag=20*log10（abs（H））;
%显示汉宁滤波器的频响
plot（f,mag）;
title（'频率响应'）;
figure（9）;
%调用滤波器对信号进行滤波
ylb=fftfilt（hn,x）;
%显示滤波后的语音信号
plot（ylb）;
title（'汉宁高通滤波后的语音信号'）;
figure（10）;
%滤波后语音信号进行FFT变换
FYLB=fft（ylb）;
%求模值
FABS=abs（FYLB（1:
N/2））;
%显示滤波后信号的频谱
plot（FABS）;
title（'汉宁高通滤波后的频谱'）;
clearall;
N1=64;
fs=10;
n=0:
1:
N1-1;
f1=2;
f2=2.05;
f3=1.9;
xn=sin（2*pi*f1*n/fs）+sin（2*pi*f2*n/fs）+sin（2*pi*f3*n/fs）;%满足条件的正弦序列x（n）
Xk=fft（xn）;
AXk=abs（Xk（1:
N1/2））;%对x（n）进行64点DFT运算
figure
（1）
plot（n,xn）;%绘出64点采样正弦序列x（n）
xlabel（'k'）;ylabel（'|X（k）|'）
figure
（2）
k=（0:
N1/2-1）*fs/N1;
plot（k,AXk）;%绘出X（k）64点DFT频谱幅度
xlabel（'k'）;ylabel（'|X（k）|'）
%X（k）补零到128点DFT
N2=128;
xn=[xnzeros（1,N2-N1）];%对64点采样序列x（n）进行补零到128点
Xk=fft（xn）;
AXk=abs（Xk（1:
N2/2））;%对补零到128点的序列进行DFT运算
m=0:
1:
N2-1;
figure（3）
k=（0:
N2/2-1）*fs/N2;
plot（k,AXk）%绘出X（k）补零到128点DFT频谱幅度
xlabel（'k'）;ylabel（'|X（k）|'）
%X（k）128点DFT
n=0:
1:
N2-1;
xn=sin（2*pi*f1*n/fs）+sin（2*pi*f2*n/fs）+sin（2*pi*f3*n/fs）;%满足条件的正弦序列
Xk=fft（xn）;
AXk=abs（Xk（1:
N2/2））;%对x（n）进行128点DFT运算
figure（4）
k=（0:
N2/2-1）*fs/N2;
plot（k,AXk）;%绘出X（k）的128点DFT频谱幅度
xlabel（'k'）;ylabel（'|X（k）|'）
64点采样正弦序列x（n）X（k）64点DFT
X（k）64点DFTX（k）补零到128点DFT
升级会员 